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Erratum: The evolution of circular, non-equatorial orbits of Kerr black holes due to 

gravitational-wave emission 
[Phys. Rev. D 61, 084004 (2000)] 

Scott A. Hughes 
PACS numbers; 04.30.Db, 04.30.-w, 04.25.Nx, 95.30.Sf 
Equation (4.52) of this paper should read 

yH.CO _ I _',\l + k yH .CO 

As a consequence, the odd / contributions to the gravitational waveforms shown are off by a minus sign. However, 
all of the rates of change discussed in this paper (i5, L^, Q, f, L) are unchanged since they only depend on the square 
of this coefficient. 
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A major focus of much current research in gravitation theory is on understanding how radia- 
tion reaction drives the evolution of a binary system, particularly in the extreme mass ratio limit. 
Such research is of direct relevance to gravitational-wave sources for space-based detectors (such as 
LISA). We present here a study of the radiative evolution of circular (i.e., constant Boyer-Lindquist 
coordinate radius), non-equatorial Kerr black hole orbits. Recent theorems have shown that, at least 
in an adiabatic evolution, such orbits evolve from one circular configuration into another, changing 
only their radius and inclination angle. This constrains the system's evolution in such a way that 
the change in its Carter constant can be deduced from knowledge of gravitational wave fluxes prop- 
agating to infinity and down the black hole's horizon. Thus, in this particular case, a local radiation 
f*^ . reaction force is not needed. In accordance with post-Newtonian weak-field predictions, we find that 

C^ ' inclined orbits radiatively evolve to larger inclination angles (although the post-Newtonian predic- 

CNj , tion overestimates the rate of this evolution in the strong field by a factor < 3). We also find that 

pH ■ the gravitational waveforms emitted by these orbits are rather complicated, particularly when the 

jrt ' hole is rapidly spinning, as the radiation is influenced by many harmonics of the orbital frequencies. 
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PACS numbers: 04.30.Db, 04.30.-w, 04.25.Nx, 95.30.Sf 

I. INTRODUCTION 
A. Motivation: the relativistic two-body problem in the extreme mass ratio limit 



^\ \ An outstanding problem in general relativity is the evolution of binary systems with strong gravity and compact 

On ' bodies. Although post-Newtonian approximations very successfully describe the evolution of such binaries when the 
(\ ' bodies are widely separated pi, no approximation schemes work well when the bodies are close: in general, there is no 
^^. small parameter that can be used to develop any approximations. Numerical relativity will be needed to accurately 
.' ' evolve and understand the dynamics of compact binaries in the strong-field regime. Despite much effort and progress 
(^JT) PI , numerical relativity is still several years away from being able to solve the most interesting strong- field problems 
ILJ such as the final inspiral and merger of binary black hole and binary neutron star systems. The full relativistic 
. ^H , two-body problem is thus likely to remain unsolved into the near future. 

j^ ' There is one limit in which the two-body problem can be solved to very high accuracy with tools available now: 

H \ the limit in which the mass of one body, mi = /i, is much less than the mass of the other, a black hole with rn2 = M . 
. .. : The body of mass fi can be treated as perturbing the "background" black hole spacetime generated by the mass Af . 
The evolution of the system can then be studied with perturbation techniques such as the Teukolsky equation [g| . 

This limit is in fact of great interest since extreme mass ratio systems {e.g., a compact body of mass 1 — lOAf0 
orbiting a black hole of mass 10^~^ -^o) are among the most important candidate sources for space-based gravitational- 
wave detectors such as the Laser Interferometer Space Antenna (LISA) Gj. Recent estimates place the number of 
such extreme mass ratio inspirals, occurring at distances D ^ 1 Gpc, in the range of 1/ycar to 1/month g,^. LISA 
should be able to measure the gravitational waves emitted by these inspirals with amplitude signal-to-noise ratios 
around 10-100 §. 

During its last year of inspiral, the compact body orbits in the very strong-field region near the massive black hole's 
event horizon, typically spiraling due to gravitational-wave emission from r ^ lOM to the innermost stable circular 
orbit. (We use units with G = c — 1.) The orbit is dynamically unstable there, so the body quickly plunges into 
the hole, cutting off the signal. During this year, the system emits roughly 10^ gravitational-wave cycles, mostly in 
LISA'S most sensitive frequency band (depending upon the mass of the black hole). If it is possible to accurately track 
the waves' phase during this year (using, for instance, the technique of matched filtering with templates), it will be 
possible to make very accurate measurements of the black hole's characteristics, and perhaps to make detailed tests 
of general relativity (for example, by "mapping" the massive object's spacetime to test whether it is in fact a black 
hole or a more exotic compact object B). 



Making such accurate measurements will require a high-precision means of modeling gravitational waves from 
extreme mass ratio binaries. One can divide in two the influences which drive such systems' evolutions: the radiation 
reaction force (backreaction due to gravitational-wave emission causing the orbiting body to lose orbital energy and 
spiral in), and environmental influences (effects due to the systems' astrophysical environment). It is possible that 
environmental influences for the most important sources will be negligible. Perhaps the most important influence is 
the interaction of the orbiting body with material accreting onto the massive black hole: black holes in the target 
mass range lO^^^M© reside at the core of galaxies where they accrete gas from their environment. In the majority 
of cases, accretion occurs at rather low rates (several orders of magnitude less than the Eddington rate [Q). For 
these "normal" galaxies, much evidence [^,0 suggests that the gas accretes via an advection dominated accretion 
flow ( ADAF) . Narayan has shown that the influence of an ADAF upon an inspiraling compact body will be far less 



important than radiation reaction |10|: the timescale for ADAF drag to change an orbit's characteristics is many 
orders of magnitude longer than the radiation reaction timescale. If accretion drag remains the most important 
enviromental influence on extreme mass ratio binaries, then radiation reaction is likely to be the only important 
element needed to construct precise models of their evolution. Because in general we expect orbits to be inclined and 
rather eccentric Q, techniques must be constructed for analyzing the radiative evolution of generic Kerr orbits. 

If the mass ratio is extreme, one can use a linear approximation to study the system's evolution. Split the spacetime 
metric into a "background" black hole piece plus a perturbation: g^p — g^0"{M, a) + hap{^). The system's evolution 
is governed by the properties of the full spacetime gap, but it is convenient to regard the motion of the small body 
dx^ /dr as a geodesic of g?§" plus corrections from an instantaneous radiation reaction force /rr(t): 



:./3 

dx^' _ dx^' 
dr dr 



odcsic 



fU^)dT. (1.1) 



The force /^j^(t) (where r is proper time measured by the orbiting body) encapsulates the manner in which the 
motion of the body deviates from motion in the background spacetime. In particular, it embodies the effects of 
radiation reaction, causing the inspiral of the small body toward the black hole as the system's orbital energy and 
angular momentum are bled off by gravitational-wave emission. Detailed quantitative understanding of this force is 
needed to precisely model the evolution of extreme mass ratio binaries. 

There currently exists a well-defined prescription for calculating the radiation reaction force. Quinn and Wald |l^] 
and Mino et al. [ [L2[ have independently and by several different techniques derived a rather general expression for 
the force, depending upon a "tail" which is integrated over the past worldline of the orbiting body's motion. This 
tail reflects the fact that, because of scatter from spacetime curvature, the domain of dependence of radiation for an 
event lies inside the light cone. Various groups are currently working on implementations of the force (see, e.g., Refs. 
|0|,n for a approaches based on direct evalulation of the Quinn- Wald-Mino quantities on the worldline, and p5|-[l^] 
for an approach based on multipole decomposition). It is likely to be some time before practical implementations for 
astrophysically realistic cases will be available. 

B. Radiation reaction without radiation reaction forces 

One can parameterize Kerr orbits by their three constants of motion: energy E^ z-component of angular momentum 
Lz, and Carter constant Q = Pg + cos^ [o^(l — E"^) + csc^ 9Lf\ |l9). A given set [E, Lz,Q) can be thought of as 
a point in the space of all possible Kerr orbits. When orbits evolve due to radiation reaction, they generate a 
trajectory [E{t),Lz{t),Q{t)] through this space. One can then regard the major goal of radiation reaction research 
as understanding all physically allowed trajectories through this orbital phase space. 

Consider the limit in which the system's evolution is adiabatic: the timescale trr for the orbit's parameters to 
change is much longer than an orbital period T. In this limit, the system spends many periods near any point on 
its phase space trajectory: each point on the trajectory is itself very nearly a geodesic orbit. [If the evolution is not 
adiabatic, the system moves along this trajectory too rapidly for this to be the case: it does not remain near any set 
of constants [E, L^, Q) long enough to approximate an orbit.] One can then regard the system as slowly passing from 
one geodesic configuration to another. 

The mechanism pushing the system from one geodesic to the next is a radiation reaction force. In this limit, the 
effect of this force can be understood as a slow change to the "constants" of the orbital motion. Might it not be 
possible to indirectly infer the change in these constants and thereby deduce the radiative evolution without actually 
computing the radiation reaction force? For ex amp le, the energy carried by gravitational waves to infinity and down 
the event horizon is well understood (cf. Refs. Ig^l ^^^ lUl)- Might it be possible to deduce from the gravitational- 
wave fiux the change of all three orbital constants and thereby deduce, in the adiabatic limit, the system's radiative 
evolution? 



In general, it is not possible to deduce the rate at which all three constants so change: the energy and the angular 
momentum can be read from the gravitational-wave flux, but usually the Carter constant cannot be. This is because 
the energy and angular momentum are scalars which are linearly constructed from an orbit's momentum p'^, whereas 
the Carter constant is a scalar which is quadratically constructed from the momentum. 

Consider a particle orbiting a Kerr black hole. At some instant, the particle has momentum pg. This particle 
radiates for some time and falls into a new orbit with instantaneous momentum p^. The change in the particle's 
momentum is Sp'^ — p^ — Pg . 

The Kerr metric admits a timelike Killing vector T^, an azimuthal Killing vector <i>^, and a Killing tensor Q^;^ (see 
p2[ for an explicit representation of this tensor in Boyer-Lindquist coordinates) . The energy of the orbiting particle 
is given hy E = —Tf^p'^. Thus, the particle has energy Eb = — T^Pg before it radiates, and energy Ea = —T^p^ after 
it radiates. The change in the energy is carried away by the radiation: 

SE^Eb-Ea^ ~T^{p'^ ~ p^) = T^<5p^ . (1.2) 

If we consider a "graviton limit"|^ of the emitted radiation, then (5pgj.aviton = —5p^ is the 4-momentum carried by 
the radiation itself. The quantity SE is the energy that observers at infinity measure the radiation to carry (or, 
that it adds to the black hole's mass after it falls down the event horizon). Since 6E = — T)j<^Pgraviton depends only 
on properties of the gravitons radiated to infinity or down the horizon, one can deduce that the particle's energy is 
changed from E^ to E^ — SE ~ E^ directly out of the gravitational- wave flux. Using $^ rather than T^, we see that 
the change in the particle's (z-component of) angular momentum can also be read from the flux. 

Consider now the particle's Carter constant. It is related to the particle's instantaneous momentum via the Killing 
tensor: Q — Q^vP^p'^ ■ Thus, before radiating, the particle has Qb = Q/^i>PbPb' ^'^'^ after has Qa — Q^ivP^p'a- 
Writing Qa using Pa ^ Pb + ^p^^ , and then evaluating Qb — Qa yields 

5Q - Qb - Qa - 2Q^,pg(5p'' + Q^Jp^Sp'^ . (1.3) 

Because this depends explicitly on the local, instantaneous momentum of the particle, this quantity cannot be read 
from the radiation flux. Indeed, dividing by St, taking the limit St — > 0, and recognizing that in the limit Sp^ /St is 
the radiation reaction force /^j^, we find that the rate of change of the Carter constant is given by 

Q = 2Q^,p^/^j, . (1.4) 

In general, the radiation reaction force is needed to compute the rate at which the Carter constant changes even in 
the adiabatic limit. 

In general, "radiation reaction without radiation reaction forces" does not work. There are, however, special cases 
where it does work. Consider first equatorial orbits. The Carter constant for equatorial orbits is zero, and, since there 
is no means by which these orbits can be raised out of the equatorial plane (the Kerr metric is reflection symmetric 
about its equator), Q remains zero at all times. Thus, the system's evolution is entirely given by the two quantities 
E and L^. This is not surprising. Equatorial orbits can be described with two parameters: a radial measure p (such 
as the semi-latus rectum; see [pl| ), and an eccentricity e. Orbital evolution is given by the rates at which these 
parameters change. Connecting the evolution of these two orbital quantities to the evolution of the two "conserved" 
quantities E and L^ fully specifies the system's evolution. Thus, equatorial orbit evolution can be fixed by examining 
radiation fiux. In particular, since Schwarzschild black holes have no preferred orientation, any Schwarzschild orbit is 
equatorial and can be evolved by studying radiation flux. Detailed studies of such orbits' evolution (and the waveforms 
they produce) can be found in |2q-|3^]. 

Kerr black holes have a preferred orientation defined by their spin axis, so equatorial and non-equatorial orbits are 
quite different. Various evolutionary studies of such orbits have been done: see [ p3| for a first analysis of waveforms and 
energy fluxes from circular equatorial orbits; |p3 for a study of waveforms and energy fluxes from eccentric equatorial 
orbits; |35[| for an analysis of the stability of circular equatorial orbits; and [[7| for an examination of circular equatorial 
orbits with an emphasis on measurability by LISA. 

Circular, non-equatorial Kerr orbits can also be analyzed without radiation reaction forces. ("Circular orbit" means 
"orbit of constant Boyer-Lindquist coordinate radius" . The properties of such orbits in the limit a — AI are discussed 



'^For the purposes of this illustrative calculation, it is far simpler to regard the radiation as a stream of particles for which we 
can identify an associated 4-momentum. In the limit in which the radiation is clearly wavelike, we must be far more careful to 
make sure that we confine our anal ysis to the wa ve zone, and to average over several wavelengths. A more careful analysis in 
this limit would likely modify Eqs. (iLSt) and (|l.4[). 



at length in |36|.) Such an orbit has non-zero Carter constant, but it turns out that, in the adiabatic hmit, its evohition 
is entirely determined by the radiated energy and angular momentum. This is because, as has recently been proved 
l9[, circular orbits remain circular as they adiabatically evolve due to radiation reaction. That is, to very high 



accuracy the system evolves from one circular configuration to another. By requiring that the circularity condition 



hold at all times, we can write down a relationship Q — Q(E,Lz) [cf. Eq. (3.5)]. Thus, measurement of the fluxes 
E and L^ is enough to entirely specify the evolution of the system. This is not surprising: inclined circular systems 
can be described by two parameters, a radius r and an inclination angle l. Despite the fact that such orbits have 
three non-trivial conserved quantities, one would guess that information about the evolution of two of these quantities 
should suffice to specify the system's evolution. 

The remainder of this paper discusses the evolution of inclined, circular orbits of Kerr black holes. A first analysis 



of such evolution was given by Shibata 1 40 1 before the "circular goes to circular" theorems were proved. Noting that 
in the a = limit the Carter constant contains information about the x and y components of the angular momentum, 
Shibata argues that even in the case a ^ the evolution of Q should be driven by L^ and Ly carried in gravitational 
radiation Q]. Neglecting certain terms since their evolutionary timescales will be much longer than other terms, 
Shibata gives an expression for Carter constant evolution in terms of quantities that can be measured in the radiation 
flux. Although we have not checked this explicitly, we suspect that Shibata's prescription is adequate for describing 
how Q changes in the weak field [where the timescale separation he describes should be quite accurate, and also where 
radiation flux down the horizon is unimportant (see pi[| )], but would not work well in the strong field. 

C. The Sasaki-Nakamura-Teukolsky formalism 



The formalism used in this paper to compute waveforms and fluxes is based on the Teukolsky equation [Eq. (4.2)]. 
The Teukolsky equation describes the (linearized) evolution of perturbations to the Kerr spacetime. In particular, it 
gives the evolution of the complex Wcyl curvature scalar ip4 = — Ca/j^^n^m'^n^m (where Ca/i^s is the Weyl curvature 
tensor and n", rh^ are legs of the Newman-Penrose null tetrad; see, e.g., [p2]). All information about the radiation 
flux at infinity and down the horizon can be extracted from ip4. 

The Teukolsky equation has a source term which depends on an integral over the world line of the orbiting particle. 
In this paper, we work in the frequency domain, thus requiring a Fourier transform of this source. For this Fourier 
transform to be formally valid, we must know the full worldline from t = — oo to oo. Since, in fact, we do not yet 
know the radiative corrections, we approximate the worldline as a Kerr geodesic orbit. This approximation is valid 
in the adiabatic limit, and is in fact how the assumption of adiabatic motion mathematically manifests itself: we use 
the zeroth order, geodesic motion to find the first order radiative correction. 

For each set of constants {E,Lz,Q) there is in fact an entire family of orbits. Each member of the family is 
distinguished by the particle's initial position at the beginning of a period. In constructing the source function, we 
assume that all of these orbits are equivalent in the sense that they generate the same radiative trajectory through 
{E,Lz,Q) parameter space. This should be a valid assumption in the adiabatic limit. To see this, pick a parameter 
space trajectory with some initial conditions, and divide it into segments of length T (a single period). In the 
adiabatic limit, each segment has identical constants, but has different initial conditions and so corresponds to a 
different member of the orbital family. Over many periods, the orbiting particle averages all members of the family. 
Thus, in the adiabatic limit, the initial conditions of the trajectory make no difference: if we had chosen different 
initial conditions, the orbiting particle would have sampled all members of the family anyway. This averaging lets us 
assume that we may pick one representative member of the orbital family and obtain reliable results. 

To actually solve the Teukolsky equation, we integrate its source over a Green's function constructed from the 
equation's homogeneous solutions. Here, we run into a well-known technical difficulty: because it has a long-ranged 
potential, the asymptotic forms of these homogeneous solutions arc somewhat ill-behaved. In a numerical computation, 
we need to be able to calculate the Teukolsky solution at some field point r by integrating from an asymptotic regime 
where its behavior is simple. Because of the long-ranged potential, it is very difficult to properly set the phase of the 
asymptotic solution: as r ^ cx), it has an outgoing piece ex e*'^'' whose amplitude grows at a rate r** times an ingoing 
piece oc e"*""" [r* is the Kerr "tortoise coordinate", cf. Eq. (^^)]. The ingoing piece is completely lost in a numerical 
calculation. 



To get around this problem, we use the Sasaki-Nakamura equation, Eq. (4.19). Solutions of the Sasaki-Nakamura 
equation are related to solutions of the Teukolsky equation by a simple transformation [Eq. (4.24)]; and, since the 
Sasaki-Nakamura equation has a short-ranged potential, its asymptotic solutions are very well-behaved. We thus 
integrate the Sasaki-Nakamura equation, perform a transformation to the Teukolsky solution, use that solution to 
construct a Green's function, and then integrate the Green's function over the Teukolsky source. From this flnal 
integration, we construct the Weyl scalar V'4 and obtain all needed details about the radiation flux. 



D. Overview and organization of this paper 



Throughout this paper, an overdot denotes d/dt and a prime denotes dj dr. An overbar denotes complex conjugation. 
Unless otherwise specified, i, r, 0, and </> refer to the Boyer-Lindquist coordinates. The superscript "rad" refers to a 
quantity carried by radiation. Hence, W'^^ is the flux of energy carried by radiation. By contrast, E is the change 
rate of an orbit's energy. We assume that E + W^'^ = ^ L^ + E'^'^. 

Section ^ gives an overview of the properties of circular Kerr geodesic orbi ts. We r eview the equations governing 
these geodesies, and then review the properties of equatorial orbits in Sec. [I A. In Sec. [IB we discuss non-equatorial 
orbits, giving relations that E, Lz, Q, and the radiu s r must satisfy for the orbit to be circular. Frequencies of motion 
for non-equatorial orbits are calculated in Sec. II C. There are two important frequencies for circular, non-equatorial 
orbits: an azimuthal frequency fi^ (connected to the time for an orbiting particle to pass through its range of (/>), and 
a polar frequency fie (related to the time for an orbiting particle to pass through its range of 9) . These frequencies are 
equal when a = 0, and approach one another in the weak-field limit. In the strong field, and particularly for rapidly 
spinning black holes, they can be quite different. 

By requiring that circular orbits remain circular 
to E and L 



In Sec. [II, we begin analyzing the effects of radiation reaction, 
as the system evolves, we derive in Sec. [II A an expression relating Q and r 



We show in Sec. IIIB 



that requiring the system to evolve adiabatically strongly constrains the applicability of this analysis to astrophysical 
systems. For orbits in the strong field , the analyses of this paper can only be considered valid if the mass ratio of the 
system is quite extreme [cf. Eq. ( 3.16| )]. 

In Sec. 1^, we describe in detail the Sasaki-Nakamura-Teukolsky formalism used to compute all radiation reaction 
quantities. Although this formalism has been used by several other authors in the past, we have found that the 
literature contains several critical errors (particularly for non-zero black hol e spin ). Thus, we discuss this formalism 
in some detail (hopefully without introducing any errors of our own) . Section IV A reviews the Teukolsky equation [pj , 
discussing the asymptotic behavior of its homogeneous solutions, using those solutions to construct a Green's function, 
and then using the Green's function to construct 'i/'4 £ind extract the waveforms and ffuxes. (Extended discussion of 
the separated 6 dependence and efficient algorithms for computing that dependence are given in Appendix H.) In Sec. 



VA| we sketch the code used to 
that this code was 



IVB, we introduce the Sasaki-Nakamura transformation and equation, and show how one can obtain the homogeneous 
Teukolsky solutions from the well-behaved Sasaki-Nakamura solution. In Sec. IV C, we develop the source term of 
the Teukolsky equation, and use it in Sec. [VE to complete the description of '04 • 

Section M describes the numerical implementation used in this study. In Sec 
compute all radiation reaction quantities. We describe in Sec. VB several tests ("sanity checks" 
required to pass: in the weak field, the code was required to reproduce results that are well known from previous 
post-Newtonian analyses; and, in the strong-field Schwarzschild limit, radiation reaction quantities for inclined orbits 
were checked to see that they obeyed a simple relation to quantities computed for the equatorial plane. (In addition, 
many quantities were spot-checked against results for equatorial orbits that were kindly provided by D. Kennefick 
and L. S. Finn; they agreed in all cases.) 

Results for strong-field Ke rr orb its are given in Sec. VI. We discuss the waveforms and energy fluxes for a few 
representative orbits in Sec. VIA. We find that as the black hole's spin increases, the effect of high harmonics 
of the orbit's fundamental frequencies becomes more important. Mathematically, this is due to the fact that the 
transmissivity of the Teukolsky potential to radiation increases as the spin parameter a is dialed up. More physically, 
it is probably due to the fact that for large spins and strong-fields, the frequencies fl^ and fig become quite different, 
and so the motion of the source becomes almost aperiodic. The orbital motion, which i s the source of the radiation, 
thus requires many harmonics of these frequencies to be accurately described. In Sec. VI B we describe a radiation 



reaction sequence for strong-field orbits near a black hole with spin parameter a = 0.8M. We examine a set of many 
orbits, parameterized by the radius r and inclination angle 6. At each parameter space point (r, t), we compute the 
direction (f, i) in which radiation reaction drives the system. These directions would be the tangent vectors to the 
trajectories [r{t),i{t)] that evolving circular orbits would follow as they evolved. An interesting result is that the 
change in inclination angle I found in the strong field is rather smaller (by a factor ^ 3) than what one would predict 
based on extrapolating post-Newtonian theory. 

The spin choice a = 0.8M is a reasonable value for a black hole whose angular momentum has been buffered by 
magnetohydrodynamic extraction of spin energy, such as are described by some models of quasar engines p3,p3 . In 
general, whenever a black hole is coupled to an accretion disk by magnetic fields, one expects there to be a buffering 
torque that can significantly impact the hole's spin |^^. However, if the hole's evolution is driven by thin-disk accretion 
(particularly by photon capture in the thin-disk limit as described in Pq|), or the hole is produced by the merger of 
two black holes (such as a result of galaxy collisions), the spin can be much larger: a — 0.998Af is the prediction of 
photon-buffered thin-disk accretion. An analysis of radiation reaction on circular orbits in this large spin regime will 
be presented in a separate paper |27]| . 



Some concluding discussion is given in Sec. VI] . In particular, we note that the assumption of adiabaticity might not 
be reasonable in the real world. This assumption would not be needed if it were possible to compute the instantaneous 
radiation reaction force /j^r. This illustrates the importance of such a force, and provides further impetus to workers 
implementing the radiation reaction force. 



II. CIRCULAR, GEODESIC KERR ORBITS 

In this section, we discuss bound circular orbits of Kerr black holes, without incorporating radiation reaction. 
Kerr geodesies are governed by the following four equations |19| : 
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(2.1a) 
(2.1b) 
(2.1c) 
(2.1d) 



The quantities E', L^, and Q ("energy", "z-component of angular momentum", and "Carter constant") specify a 
family of geodesic orbits. They are co nserved al ong any orbit in this family. Particular members of the famil y are 



s. (FT^ - ^Id|), E ee 
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specified by initial conditions. In Ec 

have been divided by /i^ and Eqs. (2.1c) and (2. Id) by fj, (where /i is the mass of the orbiting particle). Because of 



Eqs. (|2.1a|) and (|2.1b| ) 



this, E, Lz, and Q as used throughout this paper are the specific energy, angular momentum and Carter constant: 
E = £;"""'^7/z, Lz = if "''V/i, andQ = Q"''"'^'//^^. 

The right-hand side of Eq. (|2.1b ) goes to zero as the orbiting particle's position approaches the turning poi nts o f 
its 6 motion, ^max and ^min- This causes problems in a numerical implementation [cf. the source term of Sec. IV C, 
which could be written as an integral over {d9/dt)^^]. To deal with this problem, we transform to a coordinate that 
is much better behaved at the turning points. First, define the variable z — cos^ 9. Equation (p.lb|) becomes 



d9_ 
d7 



± 



± 



y/z^ [a2(i _ ^2)] ^z[Q + Ll + a2(l - E^)]TQ 
{r^ + a?-z)\/\ — z 

V/3(z+ - z){z^ -z) 
(r^ + a^z)\/l — z 



(2.2) 



The plus sign corresponds to motion from ^min to ^max, and vice versa for the minus sign. We have defined /3 ~ 
a^(l — _E^), and z± are the two roots of the quadratic in the top line of Eq. ( ^.2[ ). 

Next, define the variable x via z = z_ cos^ x- The parameter x ranges from to 27r. As x varies from to 27r, 9 
goes from ^min to 0max (at X = tt) then back to ^min- Examining dz/d9 and dz/dx we see that 



dx 
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0<x<7r; 
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\] z- - z ' 
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TT < X < 27r 



(2.3) 



Combining Eqs. ( p.2| ) and (2_^), we obtain the geodesic equation for x'- 

dx \/(3iz+ - z) 



dr 



r^ + a^z 



(2.4) 



This form is perfectly well behaved over the entire orbit. Although the interest here is in circular orbits, circularity 
was not assumed at any point in this derivation. Equation (2.4) could be useful in the study of generic Kerr orbits. 



We specialize to circular orbits from this point onward. These orbits satisfy R ^ = R'. The first condition 
guarantees that dr/dr = 0, so that the radius does not change. The second states that the orbit is eternally at a 
turning point of its radial motion. In addition, the condition R" < must be met for the orbit to be stable. 

An orbit can be fixed by specifying r and Lz'. from the conditions R = Q — R' , Q and E are determined. The 
Carter constant and angular momentum determine the amount by which the orbit is inclined from the equatorial 
plane. As in [ p8|j44[| , we will use the following inclination angle: 

cosi=^=^=. (2.5) 

This inclination angle is a constant of the motion and is very easy to compute. One could define other inclination 
angles [e.g., the arctangent of {dd / dr) / (dcj) / dr) evaluated at 6* = 7r/2, the angle at which infinitely distant observers 
would see the particle cross the equatorial plane] . We shall stick with i as defined above since it is adequate for this 
analysis. 

A. Equatorial orbits 



To begin, consider Q = 0. By Eq. (2.5) these orbits have i = 0° or 180°, lying in the equatorial plane. By Eq. 
(2.1b|), dO/dr = 0, so they remain in the equatorial plane at all times. The conditions R = = R' yield the following 



formulae for E(r) and Lz{r): 



^pro^ ; "" '1^ ^ (2.6a) 

y/1 - 3i;2 + 2qv^ 

Lr=r.l^j£^. (2.6b) 

VI - 3u2 + 2qv^ 

E-^- ] ' '"' " ^"' , (2.6c) 

Vl - 3i;2 - 2qv^ 

,t 1 + 2qv^ + q^v^ 

L, ~ —rv — , . (2.6d) 

y/1 - 3v^ ~ 2qv^ 



Here, v = -JM/r and q = a/M . The superscripts "pro" and "ret" correspond to prograde and retrograde orbits, 
respectively^. 

As seen by observers at infinity, the particle in this equatorial orbit moves with azimuthal frequency 

_d(j} _ d(l)/dT _ Af 1/2 

■^ " d^ " 'dt/d^ " ^rV^±aM^/^ ' ^^-^^ 

where the upper sign refers to prograde and the lower to retrograde orbits. 

B. Non-equatorial orbits 

Non-equatorial orbits have Q ^ 0. We will use the following algorithm to ensure that we find all stable orbits at 
some particular radius r: 

1. The most stab le orb it is t he prograde equatorial orbit, so it makes a useful starting point. Pick a value of r and 
use Eqs. ( 2.6a ) and ( ^.6b| ) to calculate this orbit's energy and angular momentum. Since it is equatorial, Q — G. 



2. Decrease the angular momentum L^, holding r fixed. Solve the system of equations R — Q — R' io find Q and 



E. These equations admit simple analytic solutions for Q{r,Lz) and E{r,Lz) [Eqs. (2.8) and (|2.9|)] 



^Throughout this paper, we distinguish between prograde and retrograde orbits by the sign of Z/2 . Another common convention 
is to switch the sign of the black hole's spin a. This is not useful here since the characteristics of the orbit should smoothly 
vary from prograde to retrograde as the orbit's inclination varies from 0° to 180°. 



3. Repeatedly decrement the angular momentum and solve again for Q{r, L^) and E{r, L^) until either the angular 
momentum reaches L^°* (indicating that we have reached the least-bound retrograde equatorial orbit) or else 
R" = 0. Orbits with R" — are marginally bound; orbits with still lower values of Lz are not stable, and are 
thus not of interest. Radiation reaction causes such orbits to catastrophically plunge into the black hole. 

As mentioned above, when r and Lz are specified, the equations R = = R' admit a simple solution for E and Q: 



E[r,Lz)-- 
Qir.LzY- 



a 



^Ll{r- M) + rA'^ 



aLzM (r2 - a?) ± Ay^r^{r - MI) + a'^r{r + M) + 0^2(^2 _ 2Mr + 2r2)] 

\{a? + r'^)Eir,Lz) - aLz] r 9 9 . n9 / ^ 9n 

^ ' ^^ '' ^ - [r' + a^E{r, Lzf - 2aE{r, Lz)Lz + hi] . 



(2.8) 
(2.9) 



There are two roots for E. Only one of these roots is physical; the other typically gives an energy less than the energy 
of the most strongly bound orbit (the prograde equatorial orbit). In this paper, we shall focus exclusively on orbits 
for which the plus sign in the denominator of Eq. (2.S) is physical. The minus sign turns out to be physical only for 
strong-field orbits of very rapidly rotating holes; an analys is of such orbits will be presented in a separate paper [ETJ. 
At any rate, since we have fixed the choice of root in Eq. (2.8), circular orbits are entirely determined by choosing r 
and Lz (and checking that R" < 0). 



C. Frequencies of non-equatorial orbits 

As the particle orbits, its motions in 9 and </> are separately periodic. The periods for these two motions are generally 
different. In this section, we derive expressions for the periods Tg and T^, and show that if we perform our analysis 
in a certain frame we need only worry about Tg. The analysis here is very similar to that of Wilkins Q, but does 
not specialize to a = M . 



1. Period of 9 motion 



We first calculate the time (as seen by observers at infinity) for the particle to move from % = to x- Combining 
Eqs. (p^) and (|1), 



dt 



7 + a'^Ez 



where 



j = E 



dx ^(3iz+ - z) 
{P + a2)2 



aLAl- 



r2 -|-a2 



(2.10) 



(2.11) 



The time it takes to go from to x is then 

toix) = / dx , 

Jo Vf3[z+ - z{x')] 



V^ 



a^E 



E 



(7r/2 - X, V^T^) - F (71/2 - X, V^W^) - E (VW^) + K (v/^TT^) 



(2.12) 



On the last line, F{(f, k) is the incomplete elliptic integral of the first kind, K{k) is the complete elliptic integral of 
the first kind; E(if, k) and E{k) are respectively the incomplete and complete elliptical integrals of the second kind 
(using the no tatio n of pT]). 

Equation ( 2.12| ) only applies to the interval < x ^ ""Z^- It is straightforward to generalize to the interval 

V2 < X < tt: 



t(x)=<o(x), 0<x<V2, 

= to(7r/2) + to(x-7r/2), 7r/2 < x < tt 



(2.13) 



Similar formulae can be written down for tt < x < 27r; it turns out that they are not needed. We could add any 
constant to t{x)', this would specify a different member of the orbital family corresponding to (_E, L^, Q). 

The particle moves through one fourth of its 6 range as x varies from to 7r/2. Hence, Tg, the period of the 
particle's motion, is given by 



Te = 4to(7r/2) = -7=K (./I^/^) + ^^'eJ^ [k (yZT^) - E [^/IZJT;) 



The corresponding frequency of 6 motion is 2-k jTg 



(2.14) 



2. Period of (j) motion 

In this subsection we calculate the angle (j) accumulated as the particle moves from to x, and use it to calculate 
the azimuthal period T^. 



Applying the transformations z = cos^ 9 — z^ cos^ x to Eq. ( 2.1c ) and combining the result with Eq. (2^) yields 

d(h 1 f L, 



dx \//3{z+ - z) \l- z 



where 



5 = aEr^-l 



A 



(2.15) 



(2.16) 



Integrating Eq. (2.15) from to x gives the amount of 4> accumulated as the particle orbits: 



where 



Mx) 



(t>(x) - Mx) , < X < 7r/2 , 

= 00(^2) + 00 (X-V2), 7r/2<x<^; 



-^= {l, [n (7r/2, -z_, Vz-/z+) - n (7r/2 - x, -Z-, Vz-/z+) 



(2.17) 



K 



(V^r7i7) - F (71/2 - X, v/^7^)] } 



Here H((^, n, k) is the incomplete elliptical integral of the third kind, again using the notation of 
the particle moves through an azimuthal angle $ given by 



$-4(/.o(7r/2) = 



VP^ 



LJV V2, -Z-, yJz^jzA - SK vW 



(2.18) 
In a period Tg, 

(2.19) 



Unless $ equals 27r (which is only the case for a = 0), the periods of 9 and (j) motion will be incommensurate. This 
is potentially problematic, since it is not clear which period is the fundamental one to use for describing the orbits 
and gravitational radiation. However, as shown by Cutler, Kennefick and Poisson (Sec. IID of pl[|), we can expand 
the 6 motion in a Fourier series of fig harmonics: 



dt 






^inflet 



(2.20) 



First, integrate this up to obtain 



(j){t) = floi + ^ 6„e'' 



tfifli 



(2.21) 



where bn = ian/nflg. Next, average this over a time Tg. The sum goes to zero by periodicity, so 

_ $ $ 

flo = 5l0 = — = —ftg . 
Ig l-K 



(2.22) 



Analyze the system in a coordinate system that is rotating at angular velocity 17^, so that cj)' = (/> — Jl^i. In this 
coordinate system, the only frequency that matters in the analysis of azimuthal motion is fig (and its harmonics): 



b'it) = ybne 



inflst 



(2.23) 



Thus f2e and its harmonics are the only frequencies that are needed in order to analyze the motion of the orbiting 
particle in this coordinate system. The only timescale we need be concerned about is Tg. 



III. RADIATION REACTION: INITIAL CONSIDERATIONS 



In this section, we begin to consider radiation reaction and how it modifies the orbits from their purely geodesic 
form. The requirement that circular orbits remain circular allows us to write down a simple condition for Q and r 
given E and Lz. This analysis is presented in Sec. Ill A. 

As discussed in the Introduction, we assume that the system evolves in an adiabatic manner: the change in any 
orbital parameter q over one period Tg must be much less than q. By insisting that the evolution be adiabatic 



everywhere, we derive conditions on the mass ratio fi/M. As we show in Sec. IIIB, this leads to stringent conditions 
when the particle is near the black hole, but essentially irrelevant conditions when the particle is in the weak field. 
When discussing astrophysical gravitational-wave sources, our results can be considered relevant only when these 
conditions are met. 



A. Circular remains circular 



As discussed in Sec. IB, it has recently been shown P7H39|] that, under adiabatic radiation reaction, circular orbits 
evolve from one circular configuration to another (except possibly in the very strong-field regime, as the orbiting 
particle begins to plunge into the black hole). An orbit which is circular remains circular. 

For the orbit to remain circular at all times, the equations i? = and i?' — must hold, where 



R = rR' + 2EEr^ 
R' = rR!' + 8EEr^ + 2 



Q 



2a^EE - 2LzLz 



Q 



Q + 2{Lz~aE){Lz-aE) 
Q + 2{Lz - aE){Lz - aE) 



Mr - a 
M 



2/n 



These quantities are found by taking the time derivative of R and R' [where R is defined in Eq. (2.1a 
Notice that i? = and R' = can be combined into a single matrix equation. 



A-VI+ B ■V2^0 



(3.1) 



(3.2) 



where 



Vl 



V2 = 



E 

Q 



(3.3) 



The matrices A and B can then be explicitly written out by gathering th e pr oper terms in E, Lz, Q, and r. The 
solution for Q and f in terms of E and Lz is now transparent: solving Eq. (3.2) for V2 yields 



V2 



= -B-^ ■A-vi = -C-vi 



(3.4) 



Thus, we may write 



Cll 



E 



Cl2 



C21 /, C22 r 



(3.5) 



where 
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cii(Q, E, L^,r) = ~AE{1 - E^)Mr^ + UEM'^r^ - 2E[a'^{l - E^) + 3{LI + Q)]Mr'^ 

+8[a^E{2 - E^) + E{Ll + Q) - 2aL,]M^r^ 

-2a[ai;[6Af2 + LI + Q + a^{l - E"^)] - 6M^L^]Mr^ 

+4a'^E[Q + {L, - aEf]M^r - Aa{L;, - aE)[Q + {L^ - aEf]M^ , (3.6a) 

ci2(g, E, L„r) = -4L,{1 - E^)Mr^ + 16(1 - E^){L, - aE)M\^ 

+2 [L,[a^{l - E^) + LI + Q]- 6M^{L, - aE)] Mr^ 

~AL, [Q + [L^ - aEf] M^r + A{L^ ~ aE) [Q + [L^ - aEf] M^ , (3.6b) 

C2i(Q, E, L^,r) = 2Er^ - 6EMr* + ia^Er^ + 2a{L, - 2aE)Mr^ + 2a-^Er - 2a^{L, - aE)M , (3.6c) 

C22(Q, E, L^,r) = 2aEMr'^ - 2a'^L^r + 2a^{L, - aE)M , (3.6d) 

d{Q, E, L,,r) = -2(1 - E'^)Mr^ + 8(1 - E^)M^r^ +[Q + LI~ 50^(1 - E"^) - eM^] Mr^ 

+2 [0^(3 - E"^) + 2aEL:, - {LI + Q)] M^r 

+2{LI + Q)M^ - iaEL.Al'^ + a^{2E^M^ ^ LI - Q)M - a^(l - E^)M . (3.6e) 

By determining E and Lz, we determine Q and r, fully fixing the evolution of the particle's orbit. In particular, the 
rate of change of the inclination angle is 

. _ d{cosL)/dt 

1 — COS"^ i 



where 

d(cos b) 



dt ^Ll + Q 






(3.8) 



B. Adiabaticity 

In this section, we impose adiabaticity on the orbiting particle and show that it constrains the mass ratio of the 
system. We first impose adiabaticity on motion in the weak field of the black hole (r ^ M), and then impose 
adiabaticity on nrotion in the strong field (r ~ M). 

1. Weak-field radiation reaction 

In the weak-field, a post-Newtonian expansion suffices to calculate the radiated energy and angular momentum. 
Ryan [p8| shows that such a post-Newtonian expansion leads to the following results for the change in the orbital 
radius and inclination angle: 

_ 64 ^ ( M^ ^ 
rwcak--y^(^— ^ 

244 PL a /Af\"/^ . 2U PL a f mV^'^ ,^ „, 

Wk = ^^^(-) sm.^— ^-(-1 .. (3.9) 



In the weak-field, Tg ~ 2'K^r'^ jM + 0{a) Q. Putting all of this together, we find to leading order 



5 M \r ^ 

4 



iwoakJe 4887r p a ( M 



15 M M \ r , 

(3.10) 



Now impose adiabaticity: the condition rwcakle/r ^ 1 leads to the condition 
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U t) 

— < 

M 1287r 



(s) 



5/2 



and Lwc&kTg/L <C 1 leads to 



fi 15 / r \ 

M "^ 4887r \m) 



(3.11) 



(3.12) 



Since, by definition, ji/M < 1/4, Eqs. (3.11) and ( 3.12| ) are always satisfied in the weak-field. Weak-field radiation 
reaction is always adiabatic, regardless of the mass ratio. This is not surprising. 



2. Strong-field radiation reaction 



In the strong- field, we do not know a priori E and L^- Nonetheless, we would like to understand what constraints 
adiabaticity places on the mass ratio, so we must at least estimate E and L^. Ryan [ p8[ gives the following quadrupole- 
order formulae, valid in the weak field: 



jTiquad 


32 /^i \2 /A'ly 


1- 


73 a fM\^^^ 

COSi 

12M \r J 


5 




r quad 


32A/ / /^ \2 /My/2 
5 \m) \r ) 


61 a fhlV'^, 
cost-h— — — (1 - 
2AM \r J ^ 


- 3 cos^ t) 



(3.13) 



(These are the rates of change of the physical energy and angular momentum, not the specific energy and angular 
momentum used in the rest of this paper.) We will take the ansatz that these formulae are valid up to factors of order 
unity even in the strong field. (Detailed analysis shows that this ansatz is reasonable; cf. Tables | - [rv| .) 

Consider prograde, equatorial orbits with a = M, r = {1 + e)M: a particle orbiting just barely outside the event 
horizon of an extreme Kerr black hole. In this limit, 



Te = M 



47r 2 

TVs 



227r 2 



0{e) 



Solving Eq. (3^) using the quadrupole formulae ( 3.13 ) with a — M , r — (1 + e)Af and expanding in e gives 

. _ 584 13657 
5V3e 15^3 

Combining these results and imposing the adiabatic condition fTg/r <C 1, we obtain 

15e2 



/^ 



< 



M 2336\/27r 



+ O(e^) 



(3.14) 



(3.15) 



(3.16) 



This is a very stringent requirement. Since our calculation will assume that adiabaticity holds over the entire evolution, 
Eq. (B.16) tells us that our results in the strong field will only be astrophysically meaningful when applied to systems 
for which the mass ratio is very extreme. (Cutler, Kennefick, and Poisson showed a similar condition holds as an 
orbiting point particle approaches the innermost stable circular orbit of the Schwarzschild spacctime pl||.) 



IV. RADIATION REACTION: THE TEUKOLSKY AND SASAKI-NAKAMURA EQUATIONS 

A. The Teukolsky equation 



We use a formalism based on the Teukolsky equation to study radiation reaction. The Teukolsky equation describes 
the behavior of the Weyl curvature component '04, which encapsulates all information about the gravitational radiation 
flux at infinity and at the event horizon. Teukolsky showed 0] that the multipolar decomposition 



V'4 



1 



dLoJ2RlrnAr)-2SZ{0) 



imcj} —iut 



(4.1) 
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separates the evolution equation for tp4. The function -2Sf^{9) is a spin- weighted spheroidal harmonic; it is discussed 
in Appendix |A[ The radial function Rimuj {f) obeys the Teukolsky equation: 



The potential is 



A^l- (^^^ ) - V{r)Ri^^{T) = -7I_(r) . (4.2) 



iiT^ + Mir - M)K , , , 

V{r) = — )r '— + SiLor + X , (4.3) 



where K = (r^ + a^)a; — ma, and|j A = £im — 2amuj + a^uP' — 2. The number £/,„ is the eigenvalue of the spheroidal 
harmonic; see Appendix H for details. 

The homogeneous Teukolsky equation admits two independent solutions, Rf^^^ and Rf^^^, with the following asymp- 
totic values: 

pH _ nholo a2 -ipmi^r* r ^, r . 

^imuj — ^imu^ e J r ^ r^ 

Din 

RL.. = BZlr'e^'^^' + ^e— '^* , r ^ oo ; (4.4) 

and 

^/™. = A™.'^V-'-*, r^oo. (4.5) 

We have introduced Pmuj — uj — mu-^-, where uj-^- = a/2Mrj^ is the "angular velocity of the horizon" (the angular 
velocity with which inertial observers at the horizon are seen to rotate due to frame dragging Bsj ) , and the "tortoise 
coordinate" 

*f \ , 2Mr+ r-r+ 2Mr_ r - r_ 
r (r) = r H In — — — In — — -;— , (4.6) 



which is derived from the rule dr* /dr — (r^ -I- a?)/ IS.. (In these relations, r± — M ± \/ M"^ — a^; recall that r+ is 
the coordinate of the event horizon.) With these solutions and using the theory of Green's functions Q], the general 
solution of the Teukolsky equation can be written 

Rir^^ir) = ztMRZM + z^MRLM , (4.7) 

where 



^ImujKr) - r,,-,.n.m noo / "'^ TTtZivl 



By construction, Z,f_(r ^ r+) = Z^^Jr ^ c^) = 0. Defining Z,^,^ ^ Z,^,^(r ^ c«), Z^^ ^ Z-^(r ^ r+), the 
asymptotic radial solutions are 

RimJr ^ r+) = Z-^Bj^^'^A^e-^'"""'-* . (4.9) 

This solution is purely ingoing at the horizon and purely outgoing at infinity, which is physical ly c orrect. It is 
convenient to absorb the factors D^^ and -B/'j°|^ into Z,^^^ and Zj^^ respectively, and rewrite Eq. (4 



^For general spin fields, A = £im — 2amu! + a^iJ^ — s{s + 1); however, we have specialized to s 
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pholc fOO poo (r'\'T ('r'\ 

From these asymptotic solutions, we construct the energy and angular momentum flux due to gravitational radiation 
that goes to infinity and down the event horizon. First, note that the particle's motion is describable as a set of 
harmonics of the frequencies f2e and ll^, and define 

Wmfe = mVl^ + kVle . (4-11) 

Then, decompose the u dependence of Z^^ and Z^^^ as 

ZlLuj = 2^ Z^^5{UJ - UJmk) , 
k 

ZZ^=Y,ZZ^,5{uj-ujmk). (4.12) 

k 

The coefficients Z^^'^ fully determine the energy and angular momentum fluxes. 

1. Fluxes as r ^ oo 

As r — + cxD, 

Mr ^ oo) ^ 2 (''+ " '^'') ■ ^^-^^^ 

The energy flux in gravitational waves, from the Isaacson stress-energy tensor [EOi, is 



[dAdt ),. ^^ - 167r\[ dt ) ^[dt 



(4.14) 



where the angle brackets denote averaging over several wavelengths. Combining Eqs. (4J.), (19), ( 4.12| ), and ( 4.14 ), 
we obtain 

dt A^oo tk 4--:^^ 



2. Fluxes as r —* r+ 

The energy and angular momentum flux at the horizon can be calculated by measuring the rate at which the event 
horizon's area increases as radiation falls into it, following the prescription of pit] as described in Eq]. The result is 



'^^^V'' -V. "^I^^^^l' urn 

^ ^^r+ Imk ™fc 

The coefficient aimk is found by transforming from Kinnerley's null tetrad (used to construct ip^) to the Hawking- 
Hartle null tetrad |^^ (which is well-behaved at the Kerr event horizon); see Q] for details. It is given by 
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aimk = TT^ [^ , (,4.1 rj 



|C; 



mk 



with e = VA'/^ - a2/4Mr+, and 

IG^feP = [(A + 2)2 + 4aw™fc - ^a^uj^kl (^^ + 36mau;„fe - SGa^w^fe) 

+ (2A + 3) {96a^ujl,k " 48maw„.fe) + 144tj^jA'f2 _ a2) _ ^4 18) 

In order to calculate all of the fluxes and from them deduce the evolution of the particle's orbit, we need a method 
to calculate the coefficients Zlf-^^f. and Z^/,. This, in turn, requires us to calculate the coefficients Bj'^/,, Sjj"'^, and 
^imk- 'To do ^O' '^^ ^^6 the Sasaki-Nakamura equation. 

B. The Sasaki-Nakamura equation 

In principle, one should be able to calculate all the necessary coefficients directly from the Teukolsky equation. 
Consider in particular -B/J^j, . From Eq. ( p~^ ) , we know that we should be able to start with a purely ingoing pulse of 
radiation at the event horizon with unit amplitude; we should then be able to use the Teukolsky equation to integrate 
out very far, and read off the ratio B"^nk/^hnk (^^ ^'^^^ ^^ -B ,°"\./ -B,^\?). This approach does not work well in a 



practical numerical implementation. The reason is that by Eq. (4.4) the outgoing solution grows with a coefficient 
relative to the ingoing coefficient — it completely swamps the ingoing solution, making it impossible to read off B^'^/^ 
with any kind of accuracy. 

The fundamental reason for this difficulty is that the Teukolsky equation's potential V{r) is long ranged. A solution 
to this difficulty was given by Sasaki and Nakamura |^2|, who discovered a transformation that takes the Teukolsky 
function R(r), governed by an equation with long-ranged potential, to the Sasaki-Nakamura function X{r), governed 
by an equation with short-ranged potential. 

The Sasaki-Nakamura equation is 

dA^ - F(r)^^ - U{r)Xi^, = . (4.19) 

dr*^ ar* 

The functions F{r) and U{r) are shown explicitly in Appendix H. Like the Teukolsky equation, the Sasaki-Nakamura 
equation admits two solutions, whose asymptotic forms are 

vH —ip-mk^* 

XLk - ^?„lP(r)e-™'''^* + At,P(r)e--'-'^*, r ^ 00 ; (4.20) 



and 



^Imk — "^Imk^ + '^Irnk^ J ' ^ '^-l- i 

XZk^Hr)e'^-'''\ r^^- (4.21) 



The function 



A B C 

P{r) = 1 + — + -— + --3 (4.22) 

allows us to m ore accurately describe the behavior of X^-°^ near infinity. Inserting X^'°° into the Sasaki-Nakamura 



equation ( 1.1£ ), and taking the limit r — > cx), we read off A, B, and C: 



A^ -Ux + 2 + 2amu) , (4.23a) 

B = -- [(A + 2)2 - (A + 2)(2 - 4amw) - 4[amw + 3iMw - amuj{amuj + 2IMuj)]] , (4.23b) 

8 

C = -- Uamijj + Bi\ - 4 + 2amuj -h SiMw) + \2(Mujf - 2A\Muo - laujflX - 3 + m^ + 2amuj)\ . (4.23c) 

6 



In the limit a ^ 0, Eqs. (4.23a)-([4.23c) reduce to the results given in 



One transforms from the Sasaki-Nakamura function to the Teukolsky function with the rule 
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r,H,oo _ ]_ 



Ar \ H.oo P Ha 



\ ] ^Imk \Xlrnk,\ 



(4.24) 



where Xinik ~ -^irnk ^1 ^J^^^^^ ■ Fi'om this rule follow relations for the various coefficients: 

/I in 

BZ.k=-^. (4.25a) 

Df^.k = —, (4.25b) 

Co 

Bf:!^ = -j^, where (4.25c) 

dirriu^ = v/2A'fr+ [(8 - 2AiMuo - l6M'^uo'^)r\ 

+ (12mm - 16M + 16amAfcj + 2AiM'^uj)r+ - 40^™^ - \2iamM + SM^] . (4.25d) 

The functions a, /3, and ry and the coefficient co are given in Appendix ^. The quantity Bf^^ is taken from Ref. [p2[ . 
The only coefficient which must be calculated numerically is A^^^^.. This is quite reasonable: from Eq. (4.2C), we 



begin with a unit amplitude, purely ingoing pulse of radiation near the event horizon and integrate outward to read 
off A™. As r ^ oo, the outgoing piece of X(l^^ remains of constant amplitude, so Aj'^^j, is easy to extract. This 
prescription is sufficient to calculate the solutions to the homogeneous Teukolsky equation. 

C. The source term 

The next part of the analysis is to compute the source of the Teukolsky equation. This term is given by p3 

Ti^Jr) =AJdndt-^{B'., + B;')_^ ^f„-(0)e-*™^e-* , (4.26) 

where the functions Bi, and B2' are 

B', = -^L-i [p-'Lo (p-'p-'T^n)] + ^^-1 [p-yj+ (p-^p-^A-ir„„-0] , (4.27) 

B;' = ^ J+ [p-'p'A-'L^, {p-'-p-'T^,.)] ^ J+ [p-'J+ (p-V-T,,„.)] . (4.28) 

Here, p ~ —!/(?' — iacos9), p — —l/{r + iacosf?) (note that p and p have the opposite sign from that used in, for 
example, [ p2| ; this is the sign convention used in, for example, ||^). The differential operators J+ and Lg are 

Lg = dg + m CSC 9 — aoj sin 9 + s cot 9 . (4.29) 

The quantities Tnm Tnm and Tmfn are projections of the stress-energy tensor for the orbiting point particle onto 
the Newman- Penrose null-tetrad legs n, m: T„„ = T^^naUf^, etc. For a point particle moving in the Kerr spacetime, 

= p [ dr u^u/^ {-gf^ 5[t- t{T)] 6[r- r{T)] 5[9 - 9{t)] 5 [^ - (/.(t)] . (4.30) 



Here, x is an arbitrary spacetime point, r is proper time measured by the moving particle, z(t) is the particle's 
worldline, and (—5)^'^ = S sin0 is the factor which converts coordinate volume to proper volume {g is the determinant 
of the covariant components of the Kerr metric.) 



As discussed in the Introduction and in Sec. [H B, the true world line of the particle is given by a Kerr geodesic plus 
radiation reaction corrections. Since we do not yet know these corrections, we use the zeroth order geodesic motion 
to write down this stress-energy tensor. Performing the integral, we find 
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T°'''{r,0,(j),t) =fi 



i°'i,P 



E sin6'i 



5[r ~ ro)S [6 ~ 0(t)] S [0 - ,^(t)] 



(4.31) 



This is specialized to circular orbits of radius tq . 

Next, project this stress-energy tensor onto the Newman-Penrose null tetrad. From 



1 /A aAsin^ei 

n« = - I — ,1,0,- 



TOn 



2 V S' ' ' E 



V2 



(iasin^, 0, E, — i(r +a ) sir 



It is useful to write the projected stress-energy tensor components in the form 

Tab = ^Q^ir ~ ro)S [9 - e{t)] S [0 - ^t)] . 



(4.32) 



(4.33) 



Setting r = 0, performing the projection, and using Eqs. ( 2.1h| )-( 2.1d| ) yields 



2 



Cy„ 



flp_ 

2Ei 



[Eir^ + a^) - ah,] 



i sin0 aE — 



L^ 



e 



i sin 6 aE 



sm 



e 



(4.34) 



The function is given in Eq. ( 2.1b ). As written, its sign is ambiguous, depending on whether the particle is ascending 
[9 decreasing) or descending {6 increasing). We take Q to be positive, and define 



C±„ = ^ [Eir^ + a^) aL^Y = C„„ , 






[E{r^ + a?) - ah,] 



i sin 6 aE 



L, 



sin" 6 



±6 






2Ei 



i sin 6 aE — 



L, 



sm 



n 2 



±e 



(4.35) 



These quantities will be used later in order to break this ambiguity when computing the source function. 

Next use Eqs. (4.27), ( 4.28| ), (4.33), and (4.34) to expand Eq. (4.26), and then repeatedly integrate by parts so that 
no 9 derivatives are taken of any delta functions. The following identity |22 3q| simplifies this integration: 



h{e)Ls [g{e)] sinOde = - / gi0)L\_^ [h{d)] sinOde , 



where 



With moderate efi^ort, we find 



L\ — de — m CSC + aw sin 6 + s cot ( 



(4.36) 



(4.37) 



TirrUr) =-4 dt e'[-*-™*(*)l i ^^ C„„ S{r - vo) l\ [p-^^t (^3^) 



+ — ^-T^ 

^72~ 

A^p^ 
H —SJ+[p '^J+{p '^pCrnmS{r- ro))] 



J+ [p^'^p^'^A^'^Cnr?iS{r-ro)] iasinOip - p)S + lIs 

Cnrn5{r - ro)Ll[p^ S dr (p-^f)] 



(4.38) 



All functions of 6 are evaluated at 0{t), the particle's 9 coordinate at t. We have written 5 as shorthand for -28^,'^ [^(0] 
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Following l32], it is very convenient to write this in the form 

Thnu> {r) = dtA'^{ [AmO + ^nmO + -4^imo] <5(r - Tq) + dr ([^„r7il + ^mml] S{r - ro)) + d^ [A,-n^,2 Sir - To)] } 



From Eq. ( 4.38 ), it is not too difRcult to work out Aat 






A2 

2V2p-^Cn 



-trt, 



-t, 



L[L!2S + 2iap sin OLl^S 
— -p-pj LlS+ [-^+P + P] iasin9S{p~ p) 



AmmO — Sp p 
. 2\ 2p C-nm 



K\ r. K ^ (K 
A +2^P^ + *^Ha 



Amml — 2,b p pC^mrn I P 



L\S + lap sai9{p — p)S 

iK 

~A 



(4.39) 

(4.40a) 
(4.40b) 

(4.40c) 

(4.40d) 

(4.40e) 
(4.40f) 



Amm2 — bp P^mm ■ 

We will r efer to quanti ties A „ „n, A „^n, A ^,^,^, et c., wh ich are the various Aau written using Eq. ( 4.35 ) rather than 
Eq. (4.34). Equations (4.34), (4.38), and (4.40a)-(4.40f ) completely specify the source for circular orbits. 



D. Evaluation of the coefficients Z^^. and Z^^ 



The next step is to substitute this source into Eq. (4.10) and integrate for the coefficients Zf^^, ^hnu- Thanks to 
the delta functions, these integrals are trivial; the result is 



Imu) 



1 



2iujB]] 



k J -oo 



Co 



rfie*'"*-"'^^*'^ i?,^„(ro) K„o + A„^o + ^mmo] 






Imu 



dr 



dr'^ 



^Imk J-oo 



^'"" ^ " SW^df^B^ I '^^ e^''^*""^^*^^ \ RZ^irn) [Ann, + A„„-.o + A^r^o] 



dR^uj 



dr 



j2 PCXD 

\A --, 4- 4 - -J + ''"'^ 



dr^ 



^7n7n2 



^rnrn2 



(4.41) 



For a numerical implementation, it is not useful to leave this in the form of an integral over infinite domain. For 
simplicity, write the above integrals as 



yH,oo ^Ha 



die*["*-™'^(*)l/^'~[ro,6'(i)] 



(4.42) 



tH.' 



To bring out the harmonic structure of the source, we would like to write the integrand in the form J^^'j^iro) exp[i(a;- 
mft^ — kfle)t\. First, define 



j"'°°[ro,0{t)] = /^'°°[ro,0(i)]e™'["**-'^Wl 



(4.43) 



From the discussion in Sec. II C , we know that 



J^'°°[^o,^(i)]=E'^™r(^o)e-'="^*(') 



(4.44) 



fc=0 



18 



where 



jH,OQ/ \ ^ 

^e Jo 



Te 



fjf 7^-00[^ i^( + \^ ^ikfldt 



Substituting for /^'°°[ro, 6'(i)] in Eq. ( [4.42D gives 



yH.i 



c 



H,oo Y^ 



Then use Eq. (4.12) to find 



k=o-'-°° 



= 2t:C 



Ha 



[ro,0{t)]i 



^^g»[^-mn^-fene]tjH,oo^^^^^ 



Ha 



fc=0 



5(U} ~ UJmk)Jyn^{ro) . 



(4.45) 



(4.46) 



Z, 



H,oc 
Irnk 



H,oo Y^ jH,oc 
fc=0 



^^C"'^l^Jrnk 



(ro) 



(4.47) 



Ha 



We must next evaluate the number Jj^j^iro). By Eqs. ( |4.43D and (f4.45| ) 



■^™r(^o) = ^1 'die'[-"-*-"^Wl/^-°°[ro,0(i)] 



(4.48) 



Because the dependence of /„^°° on i is imphcit, it is useful to change the integration variable to x- T he integrand 
then picks up a factor {dx/dt)~^ . This is well-behaved over the entire domain of the integral [cf. Eq. ( |2.10| )]. Changing 
to 9, for example, would not work well since the factor {dO/dt)^^ is singular at the turning points ^max/min- 
Writing dt — dx (dx/dt)~-^ and combining all of the results of this section, gives 



r^H _ 

^Imk ~ 



icunikTeB^ 



Imk -j- 
TTCo 



^Jo \/0\z+-z(x)\ 



''Irak 



..T.Bil,^ 



^i^rnkdiniLuTeBl^^. ^ Jq y/l3[z+ ~ z{x)] 



dx 



j + a^Ez{x) ^±.K„,t(x)-™0(x)] J. 



±h, 



where 



--^Imk VO, z{xj\ — Rlmu (^o) [^miO + ^nmO + ^mm 



dR 



H,oo 

imuj 



mmOj 



dr 



L^nml + ^mml\ 



\k[ro,z{x)] 



2 pff.oo 



d'R 



irnuj 



dr'^ 



A±^ 



7nm2 



(4.49) 
(4.50) 

(4.51) 



7H 



The coefficients Z^^ obey a useful symmetry between harmonics {m,k) and {—m,—k). From Eqs. (4.49)-(4.51) 



and by inspection of the Teukolsky equation [Eq. (4.2)], it is apparent that ^;_'^_j, = ^*"^imfc°' where a is some 
phase factor (as yet undetermined). It is simple to determine this phase factor in the Schwarzschild limit. Because of 
spherical symmetry, waveforms emitted by a particle orbiting at some angle u to the equator of a Schwarzschild hole 
and observed in the equatorial plane are equivalent to those emitted by a particle orbiting in the equatorial plane and 
observed at an angle — i. Equating the expressions for the waveforms in these two cases yields e'" = (—1)'^+', so that 



jHa 



(-1) 



''Imk 



(4.52) 



This factor in fact arises from the rotation properties of the spin- weighted spherical harmonics. Because we generate 
spheroidal harmonics as a sum of spherical harmonics, Eq. (4.52) holds for Kerr black holes as well. Hence, we need 
only consider A: > 0. 



V. NUMERICAL IMPLEMENTATION AND VALIDATION 

In this section, we first describe the structure and methods of the numerical code that was used to study circular 
orbit radiation reaction. We next present results of certain "sanity checks" that were run to make sure that, in the 
proper limits, the code reproduces well-known results. 
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A. Code implementation 

The numerical code can be broken into two pieces: a set of "harmonic" routines which calculate the coefficients 
^irnk ^'^^ ■^imk (^nd thence the energy and angular momentum fluxes), and a "driver" routine which repeatedly calls 
the harmonic routines and determines when they have "converged" (as defined below) to give the fluxes of energy 
and angular momentum emitted by that particular orbit. 

The driving algorithm has the following structure: 

1. Choose the orbital r adiu s r an d an gular momentum Lz- From this, comput e th e orbital energy E and Carter 



constant Q via Eqs. ( p.8| ) and (p.9|), and also the inclination angle t via Eq. (|2.5| ). 

2. Loop on the harmonic index ?, starting with 1 = 2, until the /-convergence criterion discussed below is satisfied. 

3. For each value of I, loop on the index m e [— /, ...,/]. 

4. For each value of m, loop on the index k until the fc-convergence criterion discussed below is satisfied. Because 



of the symmetry condition Eq. (4.52), consider only fc > 0. 



5. Compute the frequency uJmk = mil^ + kVLg using Eqs. (2.14) and (p. 19) 



6. Compute the spheroidal harmonic expansion coefficients h"-^ [cf. Eq. (A2)]. 

7. Compute the coefficients ^;,„^. This is described separately below. 

8. Compute the energy and angular momentum fluxes to inflnity and down the hole via Eqs. (4.15) and (1.16|). 



9. Check the convergence of the k loop. The convergence test implemented here is to check whether E\^^j^ < 
ek maxfc Pi^% (where Pi^% = Eg^^ + E^^, and max^ E\^^ is the largest value of E]^^ over all k values). The 
value of efc is discussed in greater detail in Sec. W\. When this condition is met for three successive k values, 
the fc-loop is ended. Otherwise, increase k by one and repeat the /c-loop. 

10. If m = Z, terminate the m-loop, else increase m by one and repeat. 

11. Check the convergence of the I loop. The convergence test implemented here is to check whether Ej^'^ < 
€i max/ Ei^'^ (where EJ^'^ = J2mk ^imk^ ^'^'^ max; Ei^'^ is the largest value of E^'^'^ over all / values). We have 
typically used e; = 10 x e^. When this condition is met for three successive I values, the Z-loop is ended. 
Otherwise, increase I by one and repeat. 

12. Compute the total change in the angular momentum and energy of the particle: 

E = -Y.Er\ Lz = -Y.Llf. (5.1) 

I I 



13. Compute the total change in the particle's radius. Carter constant and inclination angle using Eqs. (3.5)"( p.7[ ). 

14. Compute the gravitational waveform: 

h+{e, ,/>, t) - Z/.X (0, 0, = E -^zLk-2SZ-'' (0)e^(™^--*) . (5.2) 

Imk "fc 

For a given orbit with some radius and inclination angle (r, t), this algorithm gives the orbit's gravitational waveform 
and the direction (r, I) in which radiation reaction drives it to a new orbit. 

The "harmonic" routines which calculate the coefficients Z^^'^ are sufficiently involved that they merit separate 
discussion. This algorithm assumes that the orbital characteristics E, L^, Q, and r are known, as are the harmonic 
indices /, m, k (and hence the frequency uimk), and the spheroidal harmonic expansion coefficients 6"". 
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Integrate the Sasaki-Nakamura equation (4.19) inwards from "infinity" to r using Bulirsch-Stoer integration (as 
implemented in the routine bsstepO driven by odeintO from |47| , modified to integrate complex functions of 
real arguments) to get X^^.. This int egrat ion is done in Boyer-Lindquist coordinates. The value at "infinity" 



is set using the asymptotic form Eq. ( 4.21 ) . It is of course impossible to actually integrate numerically from 
r = 00, so we have implemented a variant of Richardson extrapolation p7| to accurately compute the value 
of X^i,. Let r°° be the ith "big number" which will be used to approximate infinity, and let ^°^„jj. be the 
value of Xj^f^ computed using r°° to set the r ^ oo boundary condition. We set r^ ~ 30/uJmk, i-e., roughly 5 
mode wavelengths, and set r°^^ — 2x r°° . We then construct a rational function approximation to the sequence 
-^idmk ^^ ^ function of Xi — l/r°°, and use the approximation to extrapolate to the limit a; = 0. We iterate 
until the difference in successive approximations is ~ 10"''; this typically requires 5 ~ 10 iterations. We have 
found this to be far more accurate and faster than simply setting "oo = very large number" in the code. 

Integrate the Sasaki-Nakamura equation from the event horizon to r to get Xf^^., using the asymptotic form 



Eq. ( [4.2C| ) to set the value at the horizon. We use essentially the same numerical tricks and techniques as are 
used for Xf^j.. In particular, we construct a rational function approximation to a sequence of values X^i^^. 
computed using r^ = r+ + x^, where xi ^ 1G~^ and a;i+i = Xi/2, and then extrapolate to a; = 0. (This trick is 
needed since the Sasaki-Nakamura equation is not well-behaved in Boyer-Lindquist coordinates at the horizon.) 
Roughly 4-8 iterations are needed to get a result accurate to ~ 10~^. 

3. Integrate the Sasaki-Nakamura equation from r to infinity, using X^^^. as the value at r, and thereby compute 
the coefficient A^l^f,. The same numerical techniques used to calculate X'j^j, are used here to integrate the 
equation and to reach infinity. 



^H. 



4. Compute the homogeneous Teukolsky solutions Ri^ using Eq. ([4. 241) , and Bf^f. using Eq. ( [4.25a| ). 

5. Compute z"^^ using Eqs. ( |09| ) and ( |450| ). 



This algorithm is then called by the "driver" algorithm for the step at which Z,'"^ are needed. 



B. Code validation 



In order to verify that the algorithms described in Sec. V A are working reliably, we ran a series of tests to make sure 



that the code reproduces known results in the proper limits. First, we analyzed the weak- field limit and verified that 
the code reproduces the results given in ||32|] (who analyzed gravitational radiation emitted by point particles orbiting 
black holes using a post-Newtonian expansion). Second, we took the limit a — and verified that the radiation 
emitted by a particle in an inclined orbit about a Schwarzschild black hole has the correct behavior. 

Appendices G and I of [p2| give (very long and detailed) post-Newtonian expansions of the energy flux to infinity 
and down the horizon for a particle orbiting in the equatorial plane of a Kerr black hole. These formulae allow us to 
check that the fluxes at infinity and down the horizon agree with known results as a function of the orbit's radius r and 
the black hole's spin a. We have compared both fluxes for a large number of cases and found excellent agreement (to 
the degree expected) for all parameters. Figure y is a typical example, comparing for I = 3 the numerically computed 
downhole flux for a co-rotating orbit at r = 25M as a function of black hole spin a. The fits are quite good (except for 
I = 3, m — 2) even in this relatively strong-field region because in most cases the post-Newtonian expansion formulae 
of p2[ are very robust, including many powers of M/r. (There are not as many terms given for the case / = 3, to = 2, 
so the fits are not as robust in that case.) 

These tests demonstrate the code reliably produces radiation in the Kerr equatorial plane. The next check was 
whether the code reliably produces radiation for orbits out of the equatorial plane. A simple test is to examine 
inclined orbits in the Schwarzschild limit. In this limit, the total flux radiated for a given I mode must be invariant 
as the plane of the orbit is tilted (due to spherical symmetry); however, the distribution of the radiation among the 
k an d m indices changes. The nature of this change can be deduced by examining the Teukolsky equation source, 
Eq. ( 4.26| ). In the Schwarzschild Hmit, the spin-weighted spheroidal harmonic -2S^^{9) reduces to a spin-weighted 



spherical harmonic -2^m(^) (cf. Appendix A 2). The behavior of the system under rotation follows from the behavior 
of _2yim(^) under rotation. 

The rotation behavior of zero- weight spherical harmonics {i.e., the usual well- loved spherical harmonic) is well 
understood: if one rotates about the x-axis of the coordinate system (so that the (j) angles remain unchanged) by 
some angle t, then spherical harmonics in the new (primed) coordinate system are related to those in the original 
(unprimed) coordinate system via 
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oYi^ie')^ Y, Vl,^{i)^Yi^,{e) . (5.3) 

m' — — ^ 

The function 2?m'm('^) i^ ^he Wigner D-function; an explicit expression for it is given in Eqs. (4.255)"(4.256) of B^] 
(with Arfkcn's angles a and 7 set to zero, and (3 = l). It represents one element of a matrix that rotates spherical 
harmonics. Because this matrix commutes with the differential operator d which lowers the spin weights of spherical 



harmonics (cf. Appendix A 2), Eq. (5.3) applies to any spin weight. 

From this, it is a simple matter to show that the energy emitted by a particle orbiting at angle t into some set of 
harmonic indices (/, m, k) is related to that emitted into the indices (/, m) by a particle in equatorial orbit by 



Elmk{l^) 



^nc-™).™Wr ■ (5.4) 



poq I (k—-m),m 

We have rather exhaustively compared the right- and left-hand sides of Eq. (5.4). In Figure 0, we show a typical 
comparison, / = 4, to = 2, fc = 1. The numerical points agree with the analytic D-function curve to a fractional error 
~ 10~^ — 10~^; this kind of agreement was found in all cases examined. This indicates that the off-equator capabilities 
of the code should be rehable. 

In addition to these tests, we compared this code with unpublished results from D. Kennefick (whose code was used 
for the analysis of |35)) and L. S. Finn (whose code is used for the analysis of 0). These two codes generate radiation 
for equatorial Kerr orbits. In all cases, we have found very good agreement (typically agreeing to the number of digits 
available with Finn's code, and fractional error 10^^ — 10^® compared to Kennefick's code). This, plus the validation 
tests described above, give us great confidence that the results found with this code are reliable. 

VI. RESULTS 

The major results of this analysis break into two pieces: the gravitational waves (and associated energy spectra) 
produced by particular orbits, and the sequence of orbits through which the system passes as it evolves due to radiation 
reaction. We consider these results separately below. 

A. Radiation emitted at particular orbits 



These results were computed by implementing the algorithm discussed in Sec. VA, using the value tk — 10~ (so 
that ei = 10~^). The waveforms should therefore have a fractional accuracy of about 10~^ — 10^^. We discuss in 
detail two particular strong- field orbits. 

1. r = 7M, a = 0.95M, t = 62.43° 

The gravitational waveform produced by this orbit is shown in Figure 0. To achieve fractional accuracy 10^® on 
the /-loop required summing to / — 12. For each value of /, there are of course 2/ -I- 1 values of to,; and for each value 
of m, we needed to compute 4 to 20 values of k to achieve fractional accuracy 10~^ on the fc-loop. The code uses 
roughly 2800 separate harmonics for this orbit. The fundamental orbital frequencies have values Mfl^ = 0.05424, 
MQg = 0.04954. 

A notable feature of this waveform is the presence of many short timescale features {e.g., the small bump in /i+ 
near t — 300M, and the spiky features in hx between t — llOOA/ and t — 1500M). The presence of these features 
is consistent with the breadth of this orbit's energy spectrum with respect to k (cf. Fig. Q). High frequencies play 
a rather important role in determining the radiation to infinity because the Teukolsky potential [Eq. ( |4.3| )] is rather 
transmissive to high frequency modes for large a [ p2[ . Thus, many sum and difference harmonics of the fundamental 
frequencies are needed to accurately describe the motion. The low-frequency modulation of the waveform is due to 
frame-dragging induced precession of the orbital plane — Lense-Thirring precession. 

The downhole energy spectrum for this orbit is plotted in Fig. |^. The most notable feature is that the energy 
radiated "down the hole" is negative: rather than the orbit losing energy down the horizon, this indicates that the 
orbit gains energy from the hole. This seemingly bizarre phenomenon is due to superradiant scattering ||46| ]: the 
radiation extracts energy from the ergosphere of the Kerr black hole and pumps that energy into the orbit. This 
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transfers energy from the black hole's spin to the particle's orbit, slowing the inspiral. It is essentially a manifestation 
of the Penrose process. 

Calculations such as this have excited interest in the past in the possibility of "floating orbits" |Q: orbits which 
absorb exactly as much energy as they lose to infinity, stopping their secular inspiral due to radiation loss |51| . Floating 
orbits would appear to be incredibly interesting objects: as sources of gravitational radiation, they would produce 
strong-field, nearly stationary signals, offering the possibility of very high precision studies of the Kerr ergosphere and 
the strong-field spacetime of black holes. Working with D. Kennefick, and checking all results with his code [g^, we 
have examined very strong-field circular equatorial orbits of nearly maximal Kerr black holes (for which superradiance 
is strongest, and so are most likely to have floating orbits). We find that in no cases does the energy "flowing out of 
the horizon" come close to that radiated to infinity: summing over all multipoles, we find that at best E^ ~ —0.1 E°° . 
Although we could not find this result in the published literature, other authors who have written similar codes (S. 
Detweiler and L. S. Finn in particular) have come to the same conclusion |52[| . 

Table I shows summed radiation reaction quantities for this orbit. They are compared with the values computed 
from the post-Newtonian approximation, using formulae due to Ryan ||3§|| . These numbers, coupled with the numbers 
for orbits at r = 7M, a — 0.05Af (cf. Sec. VI A 2 and Table O) indicate that the post- Newtonian formulae are not 



particularly useful in the strong-field. As predicted by Ryan, our results indicate that the inclination angle increases 
— the orbit tends to evolve toward ant i- alignment with the black hole's spin. However, the quantitative details of 
the post-Newtonian predictions are rather inaccurate. In particular, they underestimate the rate at which the radius 
changes (by about a factor of two in this case) and overestimate the rate at which the inclination angle changes (in 
this case, by a factor of three). (Interestingly, this table shows that the post-Newtonian predictions for L^ and E 
aren't nearly as inaccurate as for r and L. Small errors in L^ and E tend to magnify to large errors in f and L.) 

2. r = 7M, a = 0.05M, l = 60.17° 

The gravitational waveform produced in this case is shown in Fig. ra. Fractional accuracy of 10^^ on the ^-loop here 
required summing to I = 10. As in the case a = 0.95A/, the code needed 4 to 20 values of k for fractional accuracy 
10^^ on the A:-loop. In this case, the code needed roughly 1300 harmonics. The fundamental orbital frequencies have 
values Mn^ = 0.05407, Mfig = 0.05378. 

Two features of this waveform are noteworthy, especially in comparison with the waveform in Fig. 0. First, note 
that the low-frequency modulation of the waveform is much slower. This is because the dragging of inertial frames is 
so much slower — the modulation is caused by the plane of the particle's orbit precessing about the black hole. At 
least at lowest order, this Lense-Thirring precession frequency is proportional to the black hole's spin. Second, it is 
clear that the waveform is much "simpler" for small spin: far fewer short timescale features are present (see the lowest 
panel of Fig. ^. This is reflected in the narrowness of the spectra of energy radiated to infinity (Fig. |^). When the 
black hole's spin is small, the Teukolsky potential is not particularly transmissive to high-frequency radiation; a large 
number of sum and difference harmonics aren't needed. The radiation at infinity is largely determined by radiation 
emitted at frequencies lu ~ 2f20, cj ~ Jig -|- fi^, and uj ~ 2Qg. (For small spins, fi^ ~ Qe, so these three frequencies 
are very nearly equal.) 

The downhole energy spectrum is plotted in Fig. |[ Note that it is nowhere negative: for a — 0.05M, the particle 
does not extract any orbital energy from the black hole's spin. The ergosphere for such a slowly rotating hole is small 
and unimportant, so radiation emitted toward the hole tends to be absorbed by the event horizon rather than being 
superradiantly scattered. 

Table ^ gives the summed radiation reaction quantities, again comparing versus Ryan's post-Newtonian results. 
As in the case a = 0.95M, we see that the post-Newtonian results give the qualitatively correct result (orbital radius 
shrinks and tilt increases), but are not quantitatively accurate in the strong field. Here, post-Newtonian theory 
underestimates the rate at which the orbital radius changes by a factor of roughly three, and overestimates the rate 
at which the inclination angle changes by about 50%. 

The behavior of dE/dt summed over all / as a function of w is shown in Fig. pi For both the downhole energy and 
energy radiated to infinity, the greatest radiated flux occurs at w ~ 0.1/M ~ 2U^ ~ 217^. Peaks in the spectrum occur 
near all integer multiples of fi^ and ilg. For a — 0.95Af , these peaks are fairly broad. This is because there is power 
at many sum and difference harmonics of il^ and fig. Since the fractional difference in these frequencies is about 
10%, all power spikes are distinctly separated from one another. For a — 0.05Af , these peaks are quite narrow. The 
fractional difference between ^^ and ilg is only about 0.5% in this case, so spikes for the various sum and difference 
harmonics are not distinctly separated. 



To demonstrate convergence to the post-Newtonian results. Tables IH and IV show summed radiation reaction 



quantities for orbits at r = lOOA^, spins a = 0.95Af and a = 0.05A'/, and inclination angles near 60°. The numerical 
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results are in much better agreement with the post-Newtonian results. This is not surprising, since r — lOOM is a 
relatively weak-field region of the spacetime. It is however reassuring; these tables are further evidence that the code's 
results are trustworthy. 

B. Radiation reaction sequences 

Fig. |l^ shows a section of the parameter space in the strong-field of a Kerr black hole with spin a = 0.8A/. As 
discussed in the Introduction, this is a reasonable choice for supermassive black holes whose spin has been buffered 
by magnetohydrodynamic extraction of the holes' spin energy [P3|-p5|. Each (r, i) coordinate point in this figure 
represents an orbit. The dotted line separates stable from unstable orbits: orbits to the left of the curve are unstable 
to perturbations and catastrophically plunge into the black hole. [This line is calculated by solvin g th e system 



R = R' = R" = for the constants parameterizing marginally stable orbits, where R is defined in Eq. ( 2.1a ).] 

The effect of radiation reaction on these orbits is indicated in Fig. n^ by arrows at various points in the figure. 
The tail of each arrow indicates a particular orbit; the arrow itself is proportional to the vector [(M//z)r, (M^//i)t]. 
The arrow points in the direction to which radiation reaction pushes the particle from orbit to orbit, and its length 



indicates how rapidly radiation reaction acts. As might be expected from the discussion in Sec. VIA, the arrows in 
this figure indicate that the inclination angle does not change very quickly. One can regard the arrows as representing 
tangent vectors to a radiation reaction trajectory [r(t), i(t)]. The nearly horizontal aspect of the arrows in this figure 
indicates that radiation reaction would change l by little during a physical inspiral. 

The extremely long arrow near r = 7M, t ~ 120° in Fig. |l^ highlights the effect radiation reaction has on orbits 
which are close to the maximum inclination angle for stable orbits. That particular point is at t = 119.194°; the 
maximum inclination angle at r = 7M is t = 119.670°. Because it is extremely close to the marginally stable orbit, 
a small "push" from radiation reaction has a very large effect. Orbits that are close to the stability threshold are 
quickly pushed into the hole by gravitational-wave emission. 

Figure |l^ plots t as a function of l for various black hole spins. All curves are for orbits at r = IDA/. Notice 
that, for large spin, the distri buti on's peak is pushed away from l = 90°, where first order post-Newtonian analysis 



predicts L is maximal [cf. Eq. (3.9), which predicts i ex sin t]. One can see that the post- Newtonian prediction is more 
accurate for small spin. This is not surprising, since Eq. ( |3.9| ) is based on leading order expansions in both A-f/r and 
a/M. Although not shown here, the peak moves toward 90° as r/M increases, indicating that this shift is truly a 
strong-field effect. 

VII. CONCLUSION 

The results presented in this paper give the first strong-field radiation reaction results in which the Carter constant 
evolves in a non-trivial manner. Circular, non-equatorial orbits are probably, howeve r, th e only case in which one 
can evolve the Carter constant by examining radiation fluxes — in all other cases, Eq. (p^), relating Q to E and L^, 
will not hold. More general prescriptions for evolving the Carter constant will require calculation of an instantaneous 
radiation reaction force. Because the results presented here are constrained to circular orbits, we cannot pretend 
that they are in any way generic. However, they contain some very interesting features that may carry over to more 
general — eccentric and inclined — orbits. They also should provide useful checks to future calculations that use a 
radiation reaction force. 

The calculations presented here indicates that I is relatively small. In particular, they show that the (dimensionless) 
ratio Ml/f ^ 1. If this result holds in general, it suggests an approximation in which one holds the inclination angle 
fixed and allows the radius and eccentricity to radiatively evolve [ p3[ . Such an approximation might give a first 
cut of the trajectory a system follows through the phase space (r, e,7) of allowed orbits, perhaps serving as the first 
order solution to an iterative scheme for finding such trajectories to higher accuracy. (It is worth noting, though, 
that preliminary investigations of very strong-field orbits of rapidly rotating holes indicate that the inclination angle 
changes rather more dramatically. This result will be presented in a foUowup paper p7|.) 

The effect of spin on the gravitational waveforms (cf. the waves plotted in Figs.^ and O, and their associated 
emission spectra, Figs. ^ - |5| and [7| ~ |§) is marked: large spin causes many harmonics of the fundamental frequencies 
n^ and VL0 to influence the waveform. This is reflected by the relative breadth of emission spectra. It seems likely 
that generic orbits — which will be further influenced by a third frequency, O^ — will have this property as well, and 
we thus expect to see many harmonics of all three frequencies to influence the waveform in the general case. In terms 
of gravitational-wave observations, this impact of the spin parameter could be nice: with such a marked effect, the 
black hole's spin might be measurable to high accuracy. On the other hand, because it has such a marked effect, data 
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analysis might require an unreasonably large number of waveform templates: waveforms that are very "interesting" 
(in the sense of having a rich, detailed structure) typically require many templates so that signal-to-noise is not lost in 
measurement. To understand how many templates are needed would require constructing accurate radiation reaction 
sequences through the phase space (r, t) [(r, e, i) in the general case], building the waveforms emitted along such 
sequences, and then constructing the metric on the manifold of waveform shapes as prescribed by Owen [ p4[ . If the 
number of templates is very large, hierarchical search techniques will probably be needed. These analyses are beyond 
the scope of this paper [^ . 

Finally, the two central approximations that go into the calculations here — adiabaticity and averaging over all 
members of the orbit family — could have a large effect on gravi tational waveforms, and hence on observations by 



space based detectors such as LISA. First, as noted in Eq. ( 3.16| ), the adiabatic assumption which goes into these 
calculations requires that the mass ratio be rather extreme in the strong field. Astrophysical high-mass-ratio systems 
of interest for gravitational- wave observations are likely to have n/M ~ 10"'^ — 10~^. This might not be extreme 
enough for astrophysical systems to evolve adiabatically into the strong-field regime. However, it is impossible to relax 
the adiabatic assumption without turning to an instantaneous radiation reaction force. Second, by averaging over 
all members of the orbit family, we wash away any influence of initial conditions on the radiation reaction sequence 
— the position and momentum of the orbiting particle when observations begin. This averages the characteristics 
of the "true" waveform over the time Tioturn it takes for the orbit to return (or at the very least, come very close 
to) its initial conditions. ("True" waveform means the waveform constructed by evolving along a radiation reaction 
sequence, not just the snapshots presented here at specific moments on the sequence.) For highly eccentric, inclined 
orbits, this time could turn out to be very long. If, in particular, it turns out to be longer than the time it takes 
for radiation reaction to change the system's orbital characteristics, such averaging would not be at all accurate. In 
this case, templates constructed for the analysis of LISA-type gravitational-wave observations would require not only 
information about the trajectory through the orbital phase space (r, e, l) but also about the initial conditions of all 
possible orbits. This could drastically increase the number of needed templates. 

In the end, we find that many of the questions raised here require the instantaneous radiation reaction force 
f^Yi- Although "radiation reaction without radiation reaction force" analyses such as the one in this paper provide 
much insight and valuable information about the waveforms emitted by high-mass-ratio systems, they are essentially 
limited by the approximations that go into them, and cannot answer some of the most important questions. Because 
of the eminent need, driven by future gravitational-wave observations, to understand radiation reaction to very good 
accuracy in the high- mass-ratio limit, programs to compute the radiation reaction force should be given very high 
priority. 
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APPENDIX A: CALCULATION OF SPHEROIDAL HARMONICS 
1. Spheroidal harmonics as a sum of spherical harmonics 

The separated 9 dependence of the function ^4 is governed by the equation 



1 d f . „d^S' 



sin9- 



s^lr. 



sm0de \ de 



^ ^ , TO^ -I- 2771S COS 6* + s^ , 

(au) COS — 2aL0 s cos a — 7. + cw, 



25 



SZ^O. (Al) 



sOi 



For our purposes, we care only about the case s = —2. Solutions to this equation are the spin-weighted spheroidal 
harmonics. When acu = 0, the solutions are the spin-weighted spherical harmonics; in this case, £im = l{l + !)■ 

This fact suggests that it may be useful to expand the spheroidal harmonics in spherical harmonics. This spectral 
decomposition takes the following form: 






(A2) 



Here, it is to be understood that sYimiS) denotes the spin- weighted spherical harmonics without including the (j) 
dependence: sYiZ''''\OA) = sYim{e)e'""^ . Also, /,„;„ = max(|s|, |m|). 
It is convenient at this point to adopt Dirac-style notation, so that 

sYjm{d) -^ \sjm) , sYjm{0) -^ {sjm\ , 

sYirnie) fie) sY,m{0) sixiO dO ^ {slm\f{e)\sjm) . (A3) 



Substitute Eq. (|A2| ) into Eq. (|Al|), and use the fact that the functions sYjm satisfy (|Al[) with aui — Q and Ejm — j{j+Y). 
The result, in Dirac notation, is 



Y^ b'^'^ [(aw)2 cos2 e - 2aus cose* - j{j + 1)] \sjm) = -£i„, ^ b''^"'\sjm) . 

Next, multiply the above expression by {slm\. The various inner products are simply evaluated 

2 1 21 + I 



1 



(s/ttt-I cos 9\sjm) — -Sji + 



2j + l 



(j, 2, m, 0|/, m) (j, 2, -s, 0\l, -s) = c"\n , 



{slm\ cos 0\sjm) = J ^——(j, 1, to, 0|Z, m)(j, 1, ^ 



,01/, 



i,^i ' 



{slm\sjm) = 5ji 



(A4) 



(A5) 



The numbers (j, i, m, 0\l, n) are Clebsch-Gordan coefficients. The fact that Clebsch-Gordan coefficients appear in this 
expression greatly simplifies the sums: it tells us that c™; 2 t^ ^ only for j e [Z — 2, ^ — 1, /, / -t- 1, / -I- 2], and c™; ^ 7^ 
only for j G [/ — 1, Z, Z -f- 1]. Performing the sums, we find 

6r2(«^)'cr2,M + &ri [(«^)'cri,u - 2ac.scl^i,.i] + &r [(«^)'cO,2 - 2ac.sc0,i - U^ + 1)] 

+ 6f^i [(ac.)2cr;i,,,2 - ^aujscT+.^i^,] + htUaujfcT+2,1,2 = ~£i^br ■ (A6) 



Equation (A6) can be rewritten as a matrix equation, with the numbers bf^ representing the coefficients of the 
matrix's eigenvector, and £;,„ the matrix's eigenvalue. The matrix so defined is clearly band-diagonal, which means 
that solving for the eigenvalues and eigenvectors is rather simple. To do so, we used routines from WA. 



2. Numerical calculation of spin-weighted spherical harmonics 

The above section describes how to express the spin-weighted spheroidal harmonics as a sum of spin-weighted 
spherical harmonics; there remains the task of actually calculating the spherical harmonics. A favored method for 
stably and accurately calculating spherical harmonics of spin- weight zero is with a recurrence relation E^. It is 
relatively simple to generalize such relations to stably and accurately calculate spherical harmonics of non-zero spin 
weight. Our discussion here will be relevant to calculation of spin weights 0, —1, and —2; generalization to other spin 
weights is straightforward. See | |57| for further discussion. 

To begin, we note that the operator d lowers the spin weight of a function as follows: 



(9,y,„,(0) = -(sin( 



d_ 



{sineYsYiUff) 



^-[il + -s)il-s + l)]'/\,,r)Y,je) 



(A7) 
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The plan is to apply d to o^/m and derive formulae for the _iF;^, and then to apply d to -il^/„ and derive formulae 
for the _2^;„i- Note also that there is an operator d which raises the spin weight of a function: 



dsY,„,{e) = -(sin 



d 



[smBy^sYUS) 



6 sin 61 

[{i-s){i + s + i)Yi\,^,)Y^je) 



This operator will come in handy when computing derivatives of the spheroidal harmonics in Sec. A 3 below. 
The zero- weight spherical harmonics are written 



oYim{0) = A{l,m)Pinr{cose) , 



where 



A{l,m) = 



{21 + 1) {I -my. 



47r 



{l + my. ' 

and the associated Legendre polynomial Pim (x) can be accurately computed from the recurrence relation |47[ 

1 



Plmix) 



I 



[xi2l - l)Pl-i.rn -{l + m- l)Pl-2, 



(A8) 



(A9) 



(AlO) 



(All) 



with starting values 



Pmmix) = (-l)™(2m - 1)!!(1 - a;2)W2 



(A12) 



Now, operate on oYim with d. Combining Eqs. (A7) and (A9), and using the notation x = cos9, we see that 



-lYirJe) - - 



A{l,m) 



VT^. 



dx yr 



Plmix) . 



(A13) 



We now need recurrence relations for the functions \/l — x'^ dPim/dx and Pim/Vl — x'^. (We treat these combinations 
as fu nction s in and of themselves; this avoids problems as a; ^- ±1.) These relations are easily derived from Eqs. 
( [aTII ) and ( |aT^ ): 



Plvajx) 
\/l — x- 



l- 



x{2l- 1) 



Pi-i, 



%/r 



{l + m~ 1) 



P- 



1-2' 



VT 



^ ""^""^ - (-ir(2m - 1)!!(1 - x2)(™-i)/2 , 



Vl-x 



P 



m-\-l,7n 



(x) 



Vi~^ 



x{2m + 1) 



Prr 



^JY^^ 



(AM) 



V 1 — X- 



■dP, 



Irn 



1 



yr 



dx I — m 

dP 



{21 ~ 1) (^/]~x^ Pi-i^„, + x^l-x^ ^^^7^) - (^ + "^ - 1) v/^ 



dPi-2^ 
dx 



v/T 



dx 

■ dPm+lj. 

dx 



mx(-l)"(2TO - 1)!!(1 - x2)("^^)/2 



= (2m + 1) 



vl-a^-Pmm + x\/l 



■dP„ 



dx 



Next, operate on -li^;^ with 9. Combine Eqs. (A7) and (A13) to obtain 



2Ytje) 



A{l,m) 



^{l-l)l{l + l){l + 2) 



(l-^')x:2-2m- 



dx 



dx 



m — 2m,x 
l-a;2 



Pi^{x) 



(A15) 



(A16) 



For these harmonics, we need recurrenc e rela tions for_P;m/(l — a;^), dPim/dx, and (1 — x'^)d^Pim/dx. Again, these 
are straightforwardly derived from Eqs. (All) and (A12): 
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1 



Plniix) 



(1 — x'^) I — m 



fnl 1^ ^'-l,ni n I 1^ ^' — 2,m 

^(2^ - 1)?^ 1^ - (; + m - 1) ., ,. 



(l-a;2) 



(l-x2) 



(l-a;2) 

(1-X2) 



(-l)'"(2TO-l)!!(l-a;2)(™-2)/2 ^ 



x(2m + 1) 



Prr 



(l-a;2) 



(A17) 



dP/« 



dx I — ' 



(2/ - 1) F,_i, 



da; 



- (^ + m-l)- 



dP/- 



/-2,m 



dP„ 



dx 
dx 



mx(-l)™(2m - 1)!!(1 - a;2)(™-2)/2 
dP 

tl-i rn.7 



= (2to+1) 



dx 



dx 



(A18) 



(l-x2) 



da; 2 



1 



^ — m 



(2Z-l)(2(l-a;^)^ 



+ .(l-.2),^^^-^ 



da; 2 



(/ + m-l)(l-a;2) 



d^P;-2, 
dx2 



(1 _ ^2) ^^i™^ ^ m(-l)™(2m - 1)!!(1 - 2;2)('"-2)/2 [^^[m - 2) - (1 - x^)] , 



(l-x2) 



dx2 

j2 p 

f^ -* m+l,m 

da; 2 



= (2m + 1) 



2(1 -x^)^^ 



d2p 



(A19) 



This procedure can be continued as far as one's stamina allows; we stop here since these functions are all that is 
needed for this paper. 

3. Some derivatives of the spheroidal harmonics 

In the source term 7;,„^, we encounter the terms L^S and L\l\S [where S is shorthand for -25'!^"'= (0)]. Using 
the spectral decomposition described above, evaluating these derivatives is quite straightforward. 
Consider first L\S. Using the spectral decomposition, this becomes 



45= Y. hLl-2Yk^m{0) 



E b. 



k=l„ 



de 



sin6' 



2 cot 61 



-2yfcm(^) +awsin6' ^ bk-2Ykm{0) ■ 

fc=imin 



The operator that appears in the first term on the second line is —d [cf. Eq. (A8)] for s — —2. Thus, 

oo 

LlS^aujsm0S- ^ bk[{k - l){k + 2)]^^^ ^iYkm{0) ■ 



k=l„ 



Next consider LlL^S: 



(A20) 



(A21) 



l\lIs ^ auoLlsmOS - ^ bk[ik - l){k + 2)f^L\^iYkm{0) 
d n 



86 sine 



cot f) + au sm ( 



5in6'S' 



- Y. hu[{k-l){k + 2)] 



1/2 



d m 



cot 6 + auj sin t 



lYkm{d) , 
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de 



auj sm f 



sin 6* fine's' 



fc=imin 

d m 



auj sm W < sm 



961 sin 



sin^ 9 )■ S 



- J2 bk[ik-l)Hk + l)ik + 2)f\Ykr,^i9)-aiusme ^ bk[ik - l){k + 2)f\iYkrn{9) , 

fc=imi„ fe=imin 

d m 



89 sin f 



auj sin 6^ + 2 cot 



S 



■ ^ 5fe[(A:- l)/c(/fc + l)(fc + 2)]^/^oyfcm(6') + awsine'L^S'- (aa;sin6l)2S' 



Thus, we finally end up with 



L\lIs= Y1 bk[{k-l)k{k + l){k + 2)]^''^oYkm{0) + '2.aLos\n9LlS-{aujiiin9fS 



(A22) 



(A23) 



Eqs. ( A21 ) and ( A25 ) are then used in the source term evaluation. 



APPENDIX B: FUNCTIONS THAT APPEAR IN THE SASAKI-NAKAMURA EQUATION 



The function F{r) that appears in Eq. (4.19) is given by 

drj/dr A 



F{r) 



T] r^ + a^ 



where 



and 



r]{r) = Co + ci/r + C2/r^ + c^/r^ + c^/r^ 



Co = -l2iujM + A(A + 2) - I2auj{auj - m) , 

Ci = 8ia[3aw — \{auj — m)] , 

C2 = -2'iiaM{aLu - to) + 12a^[l - 2{auj - mf] 

C3 = 2Aia^{auj - m) - 2AM a'^ , 

a = 12a^ . 

The function U{r) that appears in Eq. ( 4.19| ) is given by 



AC/i(r) 



MG/dr 



where 



G{r) 



(^ +a ) 
2(r - Af ) rA 



+ a? (r^ + a^)^ 



C/i(r) = y(r) + ^ 



d / d/?/dr\ dr/Zdr / dp/dr 



Q = ~iK{r)/3/A^ + 3idK/dr + X + 6A/r^ 
(3 ^ 2A[-iK{r) + r-M - 2A/r] . 

The functions K{r) and V{r) are from the Teukolsky potential, Eq. (^ 



(Bl) 

(B2) 



(B3) 
(B4) 



(B5) 
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TABLE I. Comparison of certain orbital quantities computed numerically and computed using post-Newtonian expansions 
for orbits with r = IM, a = 0.95M, l = 62.43°. In this strong- field region, we do not expect post-Newtonian expressions to be 
particularly valid. Post-Newtonian theory underestimates r by about 2, and overestimates i by about 3. 



Orbital quantity 


Numerical value 


Post-Newtonian value 


{M/^jifE 


-3.3885 X 10""* 


-3.2580 X 10"^ 


{M/fi')L, 


-3.3207 X 10"^ 


-3.5926 X 10"^ 


{M/y)r 


-4.6574 X 10^'' 


-2.7499 X 10"^ 


{M''/y)i 


1.2073 X 10"'' 


3.0806 X 10-* 



TABLE IL Comparison of certain orbital quantities computed numerically and computed using post-Newtonian expansions 
for orbits with r — 7M, a — 0.05M, i = 60.17°. As in the a = 0.95M case, we do not expect post-Newtonian expressions to 
work well here in the strong-field. Post-Newtonian theory underestimates f by about 3, and overestimates i by about 50%. 



Orbital quantity 


Numerical value 


Post-Newtonian value 


{M/^^yE 


-3.9524 X 10""* 


-3.7768 X 10"^ 


(M/^^)i, 


-3.6794 X 10"^ 


-3.5210 X 10-^ 


(M//i)r 


-1.0964 X 10"' 


-3.6762 X 10"^ 


(AfV/i)^ 


1.0875 X 10"-' 


1.5867 X lO""* 



TABLE in. Comparison of certain orbital quantities computed numerically and computed using post-Newtonian expansions 
for orbits with r = lOOM, a = 0.95M, t = 60.05°. In this relatively weak-field regime, we expect to see much better agreement 
with post-Newtonian results. The numbers in this table bear out this expectation. 



Orbital quantity 


Numerical value 


Post-Newtonian value 


(M/f^rE 


-6.2194 X 10^"' 


-6.3815 X 10-"^ 


{M/^^')L, 


-3.1175 X 10"^ 


-3.1995 X 10"^ 


{M/y)r 


-1.2610 X 10"" 


-1.2733 X 10"" 


(MV/i)^ 


1.2040 X 10"'" 


1.3389 X 10"'" 



TABLE IV. Comparison of certain orbital quantities computed numerically and computed using post-Newtonian expansions 
for orbits with r = lOOM, a = 0.05M, i = 60.00°. As with the orbit at r = lOOM, a = 0.95M, we find much better agreement 
with post-Newtonian results. 



Orbital quantity 


Numerical value 


Post-Newtonian value 


{M/^,fE 


-6.2372 X 10""' 


-6.3990 X 10""' 


{M/^^')L, 


-3.1191 X 10"'^ 


-3.2000 X 10"^ 


(M//i)r 


-1.2676 X 10"-' 


-1.2797 X 10"'' 


(M7/i)^ 


6.6936 X 10"'" 


7.0439 X 10"'" 
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FIG. 1. Comparison of the flux down the event horizon for orbits at r — 25M as a function of black hole spin a/M for 
I = 3 modes. Agreement between the numerical and post-Newtonian fluxes is quite good, except for m — 2; this is because 
the post-Newtonian expansions contain many terms, and usually are quite robust. The case m = 2 is an example where the 
expansion is not as robust. The interesting upturn in the down-horizon flux for m = 1 and a > 0.5M is due to superradiant 
scattering — some incoming radiation gets scattered by the black hole's ergosphere out to infinity. 
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FIG. 2. Typical result of validation test for inclined, Schwarzschild orbits. The dotted line is the modulus squared of the 
Wigner D-function for / = 4, m = 2, A: = 1 as a function of inclination angle t; the large black points are the ratio Ei-mk{i) / E^^. 
The numerical data for the fluxes agrees with the analytical formula for the D-function to within 10^^ — 10^^. (For this plot, 
the numerical fluxes are evaluated at infinity; the results are identical to within the error when examining fluxes down the 
horizon. Also, these results are invariant — within the error bounds — as a function of orbital radius; this plot is generated 
for a particle orbiting at r = 15M.) 
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FIG. 3. The gravitational waveform produced by orbits with r = 7M, l = 62.43° about a black hole with a = 0.95M. The 
observer is in the hole's equatorial plane, 6 = 90°. The distance to the source is D. Notice that there are many sharp features 
in this waveform, indicating the strong presence of relatively large harmonics of the fundamental frequencies fi^ and He ■ This 
is consistent with the rather broad emission spectra produced by this orbit (cf. Fig. tf) . The low frequency modulation is due to 
Lense-Thirring precession (i.e., the precession of the orbital plane due to dragging of inertial frames by the black hole's spin). 
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FIG. 4. The spectrum of energy for I = 2 modes radiated to infinity by orbits with r = 7M, t = 62.43° about a black hole 
with a = 0.95M. Of particular note in this case is that the distribution is r ather broad with respect to k. This is primarily 
due to the fact that for very large spin, the Teukolsky potential [cf. Eq. (4.3)] is fairly transmissive to high frequency modes. 
Notice that it is more transmissive to corotating modes (mu! > 0) than it is to counterrotating modes: corotating modes are 
more readily scattered by the hole. 
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FIG. 5. The spectrum of energy radiated down the event horizon for I — 2 modes for orbits with r = 7M, t = 62.43° 
about a black hole with a — 0.95M. The distribution is, for the most part, sharply negative (particularly for corotating modes, 
mu} > 0), indicating superradiant scattering. This is essentially a manifestation of the Penrose process — radiation extracts 
energy from the black hole's ergosphere. 
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FIG. 6. The gravitational waveform produced by orbits with r = 7M, l — 60.14° about a black hole with a = 0.05M. The 
observer is in the hole's equatorial plane, 6 = 90°. Although not much detail regarding the waveform is visible in this figure, the 
low frequency modulation is markedly slower than in the case a = 0.95 (cf. Fig. ra). This is not surprising: the Lense-Thirring 
precession frequency is much smaller in this case since the spin is so low. The lowest panel is a zoom on h+. The time shown 
is chosen so that comparison can be easily made with the waveform for a — 0.95M, Fig. W. In contrast to the a — 0.95M case, 
the waveform is quite a bit simpler, lacking the many sharp features seen when there is large spin. This is primarily because 
the Teukolsky potential (4.3) is not nearly as transmissive to high-frequency modes when a is small. 
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FIG. 7. The spectrum of energy for I = 2 modes radiated to infinity by orbits with r = 7M, l — 60.14° about a black hole 
with a = 0.05M. The distribution is fairly narrow with respect to k, particularly when compared with the distribution for 
a = 0.95M (Fig. m. This is because the Teukolsky potential is not very transmissive to high-frequency modes for small a. Note, 
though, that the distribution shows the potential is more transmissive to corotating modes (tticj > 0) than to counterrotating 
modes, just as in the case of large a. 
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FIG. 8. The spectrum of energy radiated down the event horizon for I — 2 modes for orbits with r = 7M, i — 60.14° about 
a black hole with a = 0.05M. In this case, the distribution is nowhere negative: the ergosphere for such a slowly rotating hole 
is practically irrelevant, and as a consequence we never see any superradiant scattering. 
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FIG. 9. Energy spectra dE/dt as a function of frequency u for orbits at r = 7M. The spectra on the left are for a — 0.95M, 
i = 62.43°; the spectra on the right are for a — 0.05M, t — 60.14°. The top two spectra are dEoo/dt, the energy radiated to 
infinity; the bottom two are dEn/dt, the energy down the horizon. All spectra have been summed over I. In all cases, the 
greatest amount of radiation comes out at oj ~ 0.1/M; this corresponds to oj --^ 2^2^ ~ 2f26i- For a = 0.95M, the power is 
smeared over several frequency bins near each peak. By contrast, the power is well-confined near a single bin for a — 0.05M. In 
both cases, there is significant power at several sum and difference harmonics of fl^ and Qq. For a — 0.95M, these frequencies 
are different enough (relative difference '^ 10%) that the effect of these harmonics is quite marked. For a = 0.05M, the 
frequencies are practically identical (relative difference ~ 0.5%), so the harmonics are barely distinguishible from the main 
peak. 
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FIG. 10. The effect of radiation reaction on orbits around a black hole with a — 0.8M. The dotted line is the maximum 
allowed inclination angle; orbits tilted beyond that line are dynamically unstable and rapidly plunge into the black hole. Each 
arrow is proportional to the vector [(M//i)r, (M^//i)t]. Thus, the orientation of the arrow indicates the direction in phase 
space to which radiation reaction drives the orbit; the length of the vector indicates how strongly it is so driven. In all cases, 
the vectors point inwards and upwards — radiation reaction drives circular orbits to smaller radii and larger inclination angles 
(except when t = 0° or 180°, in which case the inclination angle does not change). The rate at which the inclination angle 
changes is rather slow, especially compared to the rate at which the radius changes. Note the very long vector (indicating 
extremely rapid orbital evolution) at t ~ 120°, r — 7M. This orbit evolves so quickly because it happens to lie extremely close 
to the maximum allowed angle — it is barely dynamically stable, so a small push from radiation reaction has marked effects. 
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FIG. 11. The rate of change of the inclination angle I versus t, parameterized by black hole spin. All curves are for orbits 
at r = lOA/. For large spin, the maximum rate of change of the angle occurs for i > 90°, in contrast to the post-Newtonian 
prediction, which yields I ex sint [cf. Eq. (B.9)]. 
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